Catégorie : code

Dans cet article, nous proposons un solveur qui s’appuie sur l’estimateur d’erreur décrit dans ce billet.

L’idée est de proposer un solveur de type SIMPLE qui, a chaque pas de temps, calcule un estimateur d’erreur et raffine ou dé-raffine le maillage selon l’importance de l’erreur liée à la discrétisation.

• files

adaptativeSimpleFoam.C

EXE = $(FOAM_APPBIN)/adaptativeSimpleFoam

• options

EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(FOAM_DEV)/ressources/errorEstimation/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -fpermissive
    

EXE_LIBS = \
    -lturbulenceModels \
    -lerrorEstimation \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -ldynamicMesh \
    -lfiniteVolume \
    -lmeshTools \
    -lfvOptions \
    -lsampling

• createFields.H

/*
 * createFileds.h
 * 
 * Copyright 2018 arep <arep@debian01>
 * 
 * This program 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.
 * 
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02110-1301, USA.
 * 
 * 
 */

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    transportProperties.lookup("nu")
);


Info<< "Creating field error\n" << endl;
volScalarField error
(
    IOobject
    (
        "error",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("error", dimensionSet(0,1,-1,0,0,0,0), Foam::scalar(0))
);

#include "createPhi.H"


label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());


singlePhaseTransportModel laminarTransport(U, phi);

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

#include "createMRF.H"<br>

• UEqn.H

    // Momentum predictor

    MRF.correctBoundaryVelocity(U);

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

• pEqn.H

{
    volScalarField rAU(1.0/UEqn.A());
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
    MRF.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);

    tmp<volScalarField> rAtU(rAU);

    if (simple.consistent())
    {
        rAtU = 1.0/(1.0/rAU - UEqn.H1());
        phiHbyA +=
            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
        HbyA -= (rAU - rAtU())*fvc::grad(p);
    }

    tUEqn.clear();

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, U, phiHbyA, rAtU(), MRF);

    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();
        }
    }

    #include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    U = HbyA - rAtU()*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);
}

• eEqn.H

/*
 * eEqn
 * 
 * Copyright 2018 arep <arep@debian01>
 * 
 * This program 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.
 * 
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02110-1301, USA.
 * 
 * 
 */
            errorEstimate<vector> eEqn
            (
                resError::div(phi, U)
              - resError::laplacian(nu, U)
             ==
               -fvc::grad(p)
            );

	    error=mag(eEqn.error());

• adaptativeSimpleFoam.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <arep@debian01>
 * 
 * This program 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.
 * 
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02110-1301, USA.

Application
    adaptativeSimpleFoam

Description
    Steady-state solver for incompressible, turbulent flow with remeshing capabilities
    using the SIMPLE algorithm.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "simpleControl.H"
#include "fvOptions.H"
#include "errorEstimate.H"
#include "resError.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
    #include "createControl.H"
    #include "createFields.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"

    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //





   // Calculate absolute flux from the mapped surface velocity
   //phi = mesh.Sf() & Uf;

    while (simple.loop())
    {


        Info<< "Time = " << runTime.timeName() << nl << endl;

	mesh.update();

        // --- Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
            #include "pEqn.H"
	    #include "eEqn.H"
        }

        laminarTransport.correct();
        turbulence->correct();


        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

Dans cet article, nous proposons une version légèrement modifiée d’un estimateur d’erreur proposé par le fork foam-extend.

Cette version de l’estimateur d’erreur recodée pour être compatible avec OpenFoam v5.0 et disponible ici l’est à titre indicatif et doit encore être validé.

• files

LIB = $(FOAM_LIBBIN)/liberrorEstimation

• option

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
    -lfiniteVolume \
    -lmeshTools

• errorEstimate.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::errorEstimate

Description
    Residual error estimation

SourceFiles
    errorEstimate.C

\*---------------------------------------------------------------------------*/

#ifndef errorEstimate_H
#define errorEstimate_H

#include "volFields.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                        Class errorEstimate Declaration
\*---------------------------------------------------------------------------*/

template<class Type>
class errorEstimate
:
    public refCount
{
    // Private data

        // Reference to GeometricField<Type, fvPatchField, volMesh>
        const GeometricField<Type, fvPatchField, volMesh>& psi_;

        //- Dimension set
        dimensionSet dimensions_;

        //- Cell residual pointer
        Field<Type> residual_;

        //- Normalisation factor
        scalarField normFactor_;


    // Private Member Functions

        //- Return boundary condition types for the error field
        wordList errorBCTypes() const;

public:

    // Static data members

    ClassName("errorEstimate");


    // Constructors

        //- Construct from components
        errorEstimate
        (
            const GeometricField<Type, fvPatchField, volMesh>& psi,
            const dimensionSet& ds,
            const Field<Type>& res,
            const scalarField& norm
        );

        //- Construct as copy
        errorEstimate(const errorEstimate<Type>&);


    // Destructor

        ~errorEstimate();


    // Member Functions

        // Access

            //- Return field
            const GeometricField<Type, fvPatchField, volMesh>& psi() const
            {
                return psi_;
            }

            //- Return residual dimensions
            const dimensionSet& dimensions() const
            {
                return dimensions_;
            }

            //- Construct and return a clone
       	    tmp<errorEstimate<Type> > clone() const;
        // Raw residual (for calculus)

            Field<Type>& res()
            {
                return residual_;
            }

            const Field<Type>& res() const
            {
                return residual_;
            }


        // Error Estimate

            //- Cell residual (volume intensive)
            tmp<GeometricField<Type, fvPatchField, volMesh> > residual() const;

            //- Normalisation factor
            tmp<volScalarField> normFactor() const;

            //- Error estimate
            tmp<GeometricField<Type, fvPatchField, volMesh> > error() const;


    // Member Operators

        void operator=(const errorEstimate<Type>&);
        void operator=(const tmp<errorEstimate<Type> >&);

        void negate();

        void operator+=(const errorEstimate<Type>&);
        void operator+=(const tmp<errorEstimate<Type> >&);

        void operator-=(const errorEstimate<Type>&);
        void operator-=(const tmp<errorEstimate<Type> >&);

        void operator+=(const GeometricField<Type,fvPatchField,volMesh>&);
        void operator+=(const tmp<GeometricField<Type,fvPatchField,volMesh> >&);

        void operator-=(const GeometricField<Type,fvPatchField,volMesh>&);
        void operator-=(const tmp<GeometricField<Type,fvPatchField,volMesh> >&);

        void operator+=(const dimensioned<Type>&);
        void operator-=(const dimensioned<Type>&);

        void operator*=(const volScalarField&);
        void operator*=(const tmp<volScalarField>&);

        void operator*=(const dimensioned<scalar>&);


    // Friend Functions

    // Friend Operators
};


// * * * * * * * * * * * * * * * Global functions  * * * * * * * * * * * * * //

template<class Type>
void checkMethod
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&,
    const char*
);

template<class Type>
void checkMethod
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&,
    const char*
);

template<class Type>
void checkMethod
(
    const errorEstimate<Type>&,
    const dimensioned<Type>&,
    const char*
);


// * * * * * * * * * * * * * * * Global operators  * * * * * * * * * * * * * //

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const dimensioned<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const dimensioned<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const dimensioned<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const dimensioned<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const dimensioned<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const dimensioned<Type>&
);


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>&,
    const tmp<errorEstimate<Type> >&
);


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>&,
    const tmp<errorEstimate<Type> >&
);


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "errorEstimate.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

• errorEstimate.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "errorEstimate.H"
#include "zeroGradientFvPatchField.H"
#include "fixedValueFvPatchField.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class Type>
Foam::wordList Foam::errorEstimate<Type>::errorBCTypes() const
{
    // Make the boundary condition type list
    // Default types get over-ridden anyway
    wordList ebct
    (
        psi_.boundaryField().size(),
        zeroGradientFvPatchField<Type>::typeName
    );

    forAll (psi_.boundaryField(), patchI)
    {
        if (psi_.boundaryField()[patchI].fixesValue())
        {
            ebct[patchI] = fixedValueFvPatchField<Type>::typeName;
        }
    }

    return ebct;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

// Construct from components
template<class Type>
Foam::errorEstimate<Type>::errorEstimate
(
    const GeometricField<Type, fvPatchField, volMesh>& psi,
    const dimensionSet& ds,
    const Field<Type>& res,
    const scalarField& norm
)
:
    psi_(psi),
    dimensions_(ds),
    residual_(res),
    normFactor_(norm)
{}


// Construct as copy
template<class Type>
Foam::errorEstimate<Type>::errorEstimate(const Foam::errorEstimate<Type>& ee)
:
    refCount(),
    psi_(ee.psi_),
    dimensions_(ee.dimensions_),
    residual_(ee.residual_),
    normFactor_(ee.normFactor_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class Type>
Foam::errorEstimate<Type>::~errorEstimate()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::errorEstimate<Type>::residual() const
{
    tmp<GeometricField<Type, fvPatchField, volMesh> > tres
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "residual" + psi_.name(),
                psi_.mesh().time().timeName(),
                psi_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            psi_.dimensions()/dimTime,
            errorBCTypes()
        )
    );

    GeometricField<Type, fvPatchField, volMesh>& res = tres();

    res.internalField() = residual_;
    res.boundaryField() == pTraits<Type>::zero;

    res.correctBoundaryConditions();

    return tres;
}


template<class Type>
Foam::tmp<Foam::volScalarField> Foam::errorEstimate<Type>::normFactor() const
{
    tmp<volScalarField> tnormFactor
    (
        new volScalarField
        (
            IOobject
            (
                "normFactor" + psi_.name(),
                psi_.mesh().time().timeName(),
                psi_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            dimless/dimTime,
            errorBCTypes()
        )
    );

    volScalarField& normFactor = tnormFactor();

    normFactor.internalField() = normFactor_;
    normFactor.boundaryField() == pTraits<Type>::zero;

    normFactor.correctBoundaryConditions();

    return tnormFactor;
}

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::errorEstimate<Type>::error() const
{
    tmp<GeometricField<Type, fvPatchField, volMesh> > tresError
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "resError" + psi_.name(),
                psi_.mesh().time().timeName(),
                psi_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            psi_.dimensions(),
            errorBCTypes()
        )
    );

    GeometricField<Type, fvPatchField, volMesh> resError = tresError();

    resError.internalField() == residual_/normFactor_;
    resError.boundaryField() = pTraits<Type>::zero;	// [ 0 0 0 ] 

    resError.correctBoundaryConditions();

    return tresError;
}
// * * * * * * * * * * * * * * *     clone         * * * * * * * * * * * * * //
//- Construct and return a clone
template<class Type>
Foam::tmp<Foam::errorEstimate<Type> > Foam::errorEstimate<Type>::clone() const
        {

	tmp<errorEstimate<Type> > terror
	(
		new errorEstimate<Type>(*this)
	);

	return terror;

        } 
// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //

template<class Type>
void Foam::errorEstimate<Type>::operator=(const Foam::errorEstimate<Type>& rhs)
{
    // Check for assignment to self
    if (this == &rhs)
    {
        FatalErrorIn
        (
            "errorEstimate<Type>::operator=(const Foam::errorEstimate<Type>&)"
        )   << "Attempted assignment to self"
            << abort(FatalError);
    }

    if (&psi_ != &(rhs.psi_))
    {
        FatalErrorIn
        (
            "errorEstimate<Type>::operator=(const errorEstimate<Type>&)"
        )   << "different fields"
            << abort(FatalError);
    }

    residual_ = rhs.residual_;
    normFactor_ = rhs.normFactor_;
}


template<class Type>
void Foam::errorEstimate<Type>::operator=(const tmp<errorEstimate<Type> >& teev)
{
    operator=(teev());
    teev.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::negate()
{
    residual_.negate();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=(const errorEstimate<Type>& eev)
{
    checkMethod(*this, eev, "+=");

    dimensions_ += eev.dimensions_;

    residual_ += eev.residual_;
    normFactor_ += eev.normFactor_;
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const tmp<errorEstimate<Type> >& teev
)
{
    operator+=(teev());
    teev.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=(const errorEstimate<Type>& eev)
{
    checkMethod(*this, eev, "+=");

    dimensions_ -= eev.dimensions_;
    residual_ -= eev.residual_;
    normFactor_ += eev.normFactor_;
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=(const tmp<errorEstimate<Type> >& teev)
{
    operator-=(teev());
    teev.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(*this, su, "+=");
    residual_ -= su.internalField();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    operator+=(tsu());
    tsu.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=
(
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(*this, su, "-=");
    residual_ += su.internalField();
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    operator-=(tsu());
    tsu.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const dimensioned<Type>& su
)
{
    residual_ -= su;
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=
(
    const dimensioned<Type>& su
)
{
    residual_ += su;
}


template<class Type>
void Foam::errorEstimate<Type>::operator*=
(
    const volScalarField& vsf
)
{
    dimensions_ *= vsf.dimensions();
    residual_ *= vsf.internalField();
    normFactor_ *= vsf.internalField();
}


template<class Type>
void Foam::errorEstimate<Type>::operator*=
(
    const tmp<volScalarField>& tvsf
)
{
    operator*=(tvsf());
    tvsf.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator*=
(
    const dimensioned<scalar>& ds
)
{
    dimensions_ *= ds.dimensions();
    residual_ *= ds.value();
    normFactor_ *= ds.value();
}


// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //

template<class Type>
void Foam::checkMethod
(
    const errorEstimate<Type>& ee1,
    const errorEstimate<Type>& ee2,
    const char* op
)
{
    if (&ee1.psi() != &ee2.psi())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, "
            "const errorEstimate<Type>&)"
        )   << "incompatible fields for operation "
            << endl << "    "
            << "[" << ee1.psi().name() << "] "
            << op
            << " [" << ee2.psi().name() << "]"
            << abort(FatalError);
    }

    if (dimensionSet::debug && ee1.dimensions() != ee2.dimensions())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, "
            "const errorEstimate<Type>&)"
        )   << "incompatible dimensions for operation "
            << endl << "    "
            << "[" << ee1.psi().name() << ee1.dimensions()/dimVolume << " ] "
            << op
            << " [" << ee2.psi().name() << ee2.dimensions()/dimVolume << " ]"
            << abort(FatalError);
    }
}


template<class Type>
void Foam::checkMethod
(
    const errorEstimate<Type>& ee,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const char* op
)
{
    if (dimensionSet::debug && ee.dimensions()/dimVolume != vf.dimensions())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, "
            "const GeometricField<Type, fvPatchField, volMesh>&)"
        )   <<  "incompatible dimensions for operation "
            << endl << "    "
            << "[" << ee.psi().name() << ee.dimensions()/dimVolume << " ] "
            << op
            << " [" << vf.name() << vf.dimensions() << " ]"
            << abort(FatalError);
    }
}


template<class Type>
void Foam::checkMethod
(
    const errorEstimate<Type>& ee,
    const dimensioned<Type>& dt,
    const char* op
)
{
    if (dimensionSet::debug && ee.dimensions()/dimVolume != dt.dimensions())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, const dimensioned<Type>&)"
        )   << "incompatible dimensions for operation "
            << endl << "    "
            << "[" << ee.psi().name() << ee.dimensions()/dimVolume << " ] "
            << op
            << " [" << dt.name() << dt.dimensions() << " ]"
            << abort(FatalError);
    }
}


// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //

namespace Foam
{

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const errorEstimate<Type>& B
)
{
    checkMethod(A, B, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() += B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const errorEstimate<Type>& B
)
{
    checkMethod(tA(), B, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() += B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(A, tB(), "+");
    tmp<errorEstimate<Type> > tC(tB.ptr());
    tC() += A;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(tA(), tB(), "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() += tB();
    tB.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const errorEstimate<Type>& B
)
{
    checkMethod(A, B, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() -= B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const errorEstimate<Type>& B
)
{
    checkMethod(tA(), B, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() -= B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(A, tB(), "-");
    tmp<errorEstimate<Type> > tC(tB.ptr());
    tC() -= A;
    tC().negate();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(tA(), tB(), "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() -= tB();
    tB.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const errorEstimate<Type>& B
)
{
    checkMethod(A, B, "==");
    return (A - B);
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const errorEstimate<Type>& B
)
{
    checkMethod(tA(), B, "==");
    return (tA - B);
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(A, tB(), "==");
    return (A - tB);
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(tA(), tB(), "==");
    return (tA - tB);
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(A, tsu(), "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(tA(), tsu(), "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, tsu(), "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), tsu(), "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(A, tsu(), "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(tA(), tsu(), "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    tC().res() -= su.internalField();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, tsu(), "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), tsu(), "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const dimensioned<Type>& su
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const dimensioned<Type>& su
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const dimensioned<Type>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const dimensioned<Type>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const dimensioned<Type>& su
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const dimensioned<Type>& su
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const dimensioned<Type>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const dimensioned<Type>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(A, su, "==");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(tA(), su, "==");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(A, tsu(), "==");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(tA(), tsu(), "==");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const dimensioned<Type>& su
)
{
    checkMethod(A, su, "==");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const dimensioned<Type>& su
)
{
    checkMethod(tA(), su, "==");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField& vsf,
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() *= vsf;
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>& tvsf,
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() *= tvsf;
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField& vsf,
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() *= vsf;
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>& tvsf,
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() *= tvsf;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>& ds,
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() *= ds;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>& ds,
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() *= ds;
    return tC;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //

• resError.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Namespace
    Foam::resError

Description
    Namespace for residual error estimate operators.

\*---------------------------------------------------------------------------*/

#ifndef resError_H
#define resError_H

#include "resErrorDiv.H"
#include "resErrorLaplacian.H"
#include "resErrorSup.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

• resErrorDiv.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

InNamespace
    Foam::resError

Description
    Residual error estimate for the fv convection operators.

SourceFiles
    resErrorDiv.C

\*---------------------------------------------------------------------------*/

#ifndef resErrorDiv_H
#define resErrorDiv_H

#include "errorEstimate.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace resError
{
    // Divergence terms

        template<class Type>
        tmp<errorEstimate<Type> > div
        (
            const surfaceScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > div
        (
            const tmp<surfaceScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

}  // End namespace resError


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "resErrorDiv.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

• resErrorDiv.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description
    Residual error estimate for the fv convection operators.

\*---------------------------------------------------------------------------*/

#include "resErrorDiv.H"
#include "fvc.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace resError
{

template<class Type>
tmp<errorEstimate<Type> >
div
(
    const surfaceScalarField& flux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    const scalarField& vols = mesh.V();
    const surfaceVectorField& faceCentres = mesh.Cf();
    const volVectorField& cellCentres = mesh.C();
    const fvPatchList& patches = mesh.boundary();
    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    Field<Type> res(vols.size(), pTraits<Type>::zero);
    scalarField aNorm(vols.size(), 0.0);

    // Get sign of flux
    const surfaceScalarField signF = pos(flux);

    // Calculate gradient of the solution
    // Change of return type due to gradient cacheing.  HJ, 22/Apr/2016
    const tmp
    <
        GeometricField
        <
            typename outerProduct<vector, Type>::type,
            fvPatchField,
            volMesh
        >
    > tgradVf = fvc::grad(vf);

    const GeometricField
    <
        typename outerProduct<vector, Type>::type,
        fvPatchField,
        volMesh
    >& gradVf = tgradVf();

    // Internal faces
    forAll (owner, faceI)
    {
        // Calculate the centre of the face
        const vector& curFaceCentre = faceCentres[faceI];

        // Owner
        vector ownD = curFaceCentre - cellCentres[owner[faceI]];

        // Subtract convection
        res[owner[faceI]] -=
        (
            vf[owner[faceI]]
          + (ownD & gradVf[owner[faceI]])
        )*flux[faceI];

        aNorm[owner[faceI]] += signF[faceI]*flux[faceI];

        // Neighbour
        vector neiD = curFaceCentre - cellCentres[neighbour[faceI]];

        // Subtract convection
        res[neighbour[faceI]] +=
        (
            vf[neighbour[faceI]]
          + (neiD & gradVf[neighbour[faceI]])
        )*flux[faceI];

        aNorm[neighbour[faceI]] -= (1.0 - signF[faceI])*flux[faceI];
    }

    forAll (patches, patchI)
    {
        const vectorField& patchFaceCentres =
            faceCentres.boundaryField()[patchI];

        const scalarField& patchFlux = flux.boundaryField()[patchI];
        const scalarField& patchSignFlux = signF.boundaryField()[patchI];

        const labelList& fCells = patches[patchI].faceCells();

        forAll (fCells, faceI)
        {
            vector d =
                patchFaceCentres[faceI] - cellCentres[fCells[faceI]];

            // Subtract convection
            res[fCells[faceI]] -=
            (
                vf[fCells[faceI]]
              + (d & gradVf[fCells[faceI]])
            )*patchFlux[faceI];

            aNorm[fCells[faceI]] += patchSignFlux[faceI]*patchFlux[faceI];
        }
    }

    res /= vols;
    aNorm /= vols;

    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            flux.dimensions()*vf.dimensions(),
            res,
            aNorm
        )
    );
}


template<class Type>
tmp<errorEstimate<Type> >
div
(
    const tmp<surfaceScalarField>& tflux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Div(resError::div(tflux(), vf));
    tflux.clear();
    return Div;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace resError

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

• resErrorLaplacian.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

InNamespace
    Foam::resError

Description
    Residual error estimate for the fv laplacian operators

SourceFiles
    resErrorLaplacian.C

\*---------------------------------------------------------------------------*/

#ifndef resErrorLaplacian_H
#define resErrorLaplacian_H

#include "errorEstimate.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace resError
{
    // Laplacian terms

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const dimensionedScalar&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const volScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<volScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const surfaceScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<surfaceScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const volTensorField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<volTensorField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const surfaceTensorField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<surfaceTensorField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "resErrorLaplacian.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

• resErrorLaplacian.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description
    Residual error estimate for the fv laplacian operators.

\*---------------------------------------------------------------------------*/

#include "resErrorLaplacian.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace resError
{

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    surfaceScalarField Gamma
    (
        IOobject
        (
            "gamma",
            vf.time().constant(),
            vf.db(),
            IOobject::NO_READ
        ),
        vf.mesh(),
        dimensionedScalar("1", dimless, 1.0)
    );

    return resError::laplacian(Gamma, vf);
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const dimensionedScalar& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    surfaceScalarField Gamma
    (
        IOobject
        (
            gamma.name(),
            vf.time().timeName(),
            vf.db(),
            IOobject::NO_READ
        ),
        vf.mesh(),
        gamma
    );

    return resError::laplacian(Gamma, vf);
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const volScalarField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return resError::laplacian(fvc::interpolate(gamma), vf);
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<volScalarField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Laplacian(resError::laplacian(tgamma(), vf));
    tgamma.clear();
    return Laplacian;
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const surfaceScalarField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    const scalarField& vols = mesh.V();
    const surfaceVectorField& Sf = mesh.Sf();
    const surfaceScalarField magSf = mesh.magSf();
    const fvPatchList& patches = mesh.boundary();
    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    const surfaceScalarField& delta =
        mesh.surfaceInterpolation::deltaCoeffs();

    Field<Type> res(vols.size(), pTraits<Type>::zero);
    scalarField aNorm(vols.size(), 0.0);

    // Calculate gradient of the solution.
    // Change of return type due to gradient cacheing.  HJ, 22/Apr/2016
    const tmp
    <
        GeometricField
        <
            typename outerProduct<vector, Type>::type,
            fvPatchField,
            volMesh
        >
    > tgradVf = fvc::grad(vf);

    const GeometricField
    <
        typename outerProduct<vector, Type>::type,
        fvPatchField,
        volMesh
    >& gradVf = tgradVf();

    // Internal faces
    forAll (owner, faceI)
    {
        // Owner

        // Subtract diffusion
        res[owner[faceI]] -=
            gamma[faceI]*(Sf[faceI] & gradVf[owner[faceI]]);

        aNorm[owner[faceI]] += delta[faceI]*gamma[faceI]*magSf[faceI];

        // Neighbour

        // Subtract diffusion
        res[neighbour[faceI]] +=
            gamma[faceI]*(Sf[faceI] & gradVf[neighbour[faceI]]);

        aNorm[neighbour[faceI]] += delta[faceI]*gamma[faceI]*magSf[faceI];

    }

    forAll (patches, patchI)
    {
        const vectorField& patchSf = Sf.boundaryField()[patchI];
        const scalarField& patchMagSf = magSf.boundaryField()[patchI];
        const scalarField& patchGamma = gamma.boundaryField()[patchI];
        const scalarField& patchDelta = delta.boundaryField()[patchI];

        const labelList& fCells = patches[patchI].faceCells();

        forAll (fCells, faceI)
        {
            // Subtract diffusion
            res[fCells[faceI]] -=
                patchGamma[faceI]*
                (
                    patchSf[faceI] & gradVf[fCells[faceI]]
                );

            aNorm[fCells[faceI]] +=
                patchDelta[faceI]*patchGamma[faceI]*patchMagSf[faceI];
        }
    }

    res /= vols;
    aNorm /= vols;

    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            delta.dimensions()*gamma.dimensions()*magSf.dimensions()
            *vf.dimensions(),
            res,
            aNorm
        )
    );
}

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<surfaceScalarField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > tresError(resError::laplacian(tgamma(), vf));
    tgamma.clear();
    return tresError;
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const volTensorField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    return resError::laplacian
    (
        (mesh.Sf() & fvc::interpolate(gamma) & mesh.Sf())
        /sqr(mesh.magSf()),
        vf
    );
}

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<volTensorField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Laplacian = resError::laplacian(tgamma(), vf);
    tgamma.clear();
    return Laplacian;
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const surfaceTensorField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    return resError::laplacian
    (
        (mesh.Sf() & gamma & mesh.Sf())/sqr(mesh.magSf()),
        vf
    );
}

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<surfaceTensorField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Laplacian = resError::laplacian(tgamma(), vf);
    tgamma.clear();
    return Laplacian;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace resError

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

• resErrorSup.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

InNamespace
    Foam::resError

Description
    Residual error estimate for the fv source operators

SourceFiles
    resErrorSup.C

\*---------------------------------------------------------------------------*/

#ifndef resErrorSup_H
#define resErrorSup_H

#include "errorEstimate.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{


namespace resError
{
    // Implicit source

        template<class Type>
        tmp<errorEstimate<Type> > Sp
        (
            const volScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > Sp
        (
            const tmp<volScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );


        template<class Type>
        tmp<errorEstimate<Type> > Sp
        (
            const dimensionedScalar&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );


    // Implicit/Explicit source depending on sign of coefficient

        template<class Type>
        tmp<errorEstimate<Type> > SuSp
        (
            const volScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > SuSp
        (
            const tmp<volScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "resErrorSup.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

• resErrorSup.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description
    Residual error estimate for the fv source operators.

\*---------------------------------------------------------------------------*/

#include "resErrorSup.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace resError
{

template<class Type>
tmp<errorEstimate<Type> >
Sp
(
    const volScalarField& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            sp.dimensions()*vf.dimensions()*dimVolume,
            sp.internalField()*vf.internalField(),
            scalarField(vf.internalField().size(), 0)
        )
    );
}

template<class Type>
tmp<errorEstimate<Type> >
Sp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > tee = resError::Sp(tsp(), vf);
    tsp.clear();
    return tee;
}


template<class Type>
tmp<errorEstimate<Type> >
Sp
(
    const dimensionedScalar& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            sp.dimensions()*vf.dimensions()*dimVolume,
            sp.value()*vf.internalField(),
            scalarField(vf.internalField().size(), 0)
        )
    );
}


template<class Type>
tmp<errorEstimate<Type> >
SuSp
(
    const volScalarField& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return Sp(sp, vf);
}

template<class Type>
tmp<errorEstimate<Type> >
SuSp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > tee = resError::SuSp(tsp(), vf);
    tsp.clear();
    return tee;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace resError

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //