Auteur : Alexis Sauvageon

code

OpenFOAM – adaptativeSimpleFoam

Alexis Sauvageon

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;
}


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

 

 

code

OpenFOAM – errorEstimate

Alexis Sauvageon

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

Il s’agit d’un estimateur d’erreur en résidu qui permet d’estimer l’erreur de discrétisation due à la méthode des volumes finis
sur les éléments d’un maillage 2D ou 3D. Pour plus d’informations sur les estimateurs d’erreur, on pourra se reporter aux travaux de H. Jsak ou encore dans le secteur de la mécanique des milieux continus à la documentation de Code_Aster.

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

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


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Non classé

OpenFOAM – fixedTabulatedVectorValue

Alexis Sauvageon

 

Cette condition limite permet d’imposer au bord du domaine fluide des valeurs vectorielles stockées dans une table. Les valeurs vectorielles doivent être associées à un point définit par ses coordonnées.

xyzUxUyUz
-12-420 -0.0124975-0.4567420.00380209
-12-420.274078-0.0248316-0.7279960.00482216
-12-41.60640.274078-0.011896 -0.456030.00395441
⋮⋮⋮⋮⋮⋮

Le programme utilise un algorithme des plus proches voisins pour associer à chaque centre de faces du patch sur lequel est imposée la condition limite une valeur d’entrée calculée par la méthode de la pondération inverse à la distance.

Les développements font appel à la librarie nanoflann accessible en cliquant sur le lien suivant. Cette librarie ne requiert aucune compilation préalable, toutefois les liens vers les headers doivent être correctement définis pour pouvoir compiler les sources de fixedTabularVectorValue.

Dans la suite, nous fournissons les sources du .h et du .cpp ainsi que des deux fichiers options et files nécessaires à la compilation. Les chemins sont susceptibles de changer en fonction de l’environnement de travail.

 

  • options
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(FOAM_DEV)/libraries/nanoflann/include \
-I$(FOAM_DEV)/libraries/utils/include
EXE_LIBS =<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • files
fvPatchFields = ./
derivedFvPatchFields = $(fvPatchFields)/derived

$(derivedFvPatchFields)/fixedTabulatedVectorValue/fixedTabulatedVectorValueFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/libArepCustomPatch<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • fixedTabulatedVectorValueFvPatchScalarField.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Class
    fixedTabulatedVectorValueFvPatchVectorField

Description
    Compute boundaries conditions from a field stored in table (format .raw)
	[X Y Z VX VY VZ]
    the two first lines of the table are not interpreted.
    The interpollated value is computed thanks to a PID method.

    \heading Patch usage

    \table
        Property    | Description             		| Required   	| Default value
        fileDir     | path to the input file  		| yes        	| 
        nbPoint     | number of point for PID 	    | no         	| 4
        wFactor     | p factor for PID      		| no        	| 2.0

    \endtable

    Example of the boundary condition specification:
    \verbatim
    myPatch
    {
        type            fixedTabulatedVectorValue;
        fileDir         "constant/inputTableData/p_boxValue.raw";
        nbPoint              8;
        wFactor            2.0;

    }
    \endverbatim


SourceFiles
    fixedTabulatedVectorValueFvPatchVectorField.C

Author
    Alexis Sauvageon.  All rights reserved

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

#ifndef fixedTabulatedVectorValueFvPatchVectorField_H
#define fixedTabulatedVectorValueFvPatchVectorField_H

#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"

#include <fstream>
#include <sstream>
#include <iterator>



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

namespace Foam
{

/*---------------------------------------------------------------------------*\
              Class fixedTabulatedVectorValueFvPatchField Declaration
\*---------------------------------------------------------------------------*/

class fixedTabulatedVectorValueFvPatchVectorField
:
    public fixedValueFvPatchVectorField
{
    // Private data

        //- Peak velocity magnitude
        std::string fileDir_;
	unsigned int nbPoint_;
	scalar wFactor_;

public:

    //- Runtime type information
    TypeName("fixedTabulatedVectorValue");


    // Constructors

        //- Construct from patch and internal field
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fvPatch&,
            const DimensionedField<vector, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fvPatch&,
            const DimensionedField<vector, volMesh>&,
            const dictionary&
        );

        //- Construct by mapping given fixedTabulatedVectorValueFvPatchVectorField
        //  onto a new patch
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fixedTabulatedVectorValueFvPatchVectorField&,
            const fvPatch&,
            const DimensionedField<vector, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchVectorField> clone() const
        {
            return tmp<fvPatchVectorField>
            (
                new fixedTabulatedVectorValueFvPatchVectorField(*this)
            );
        }

        //- Construct as copy setting internal field reference
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fixedTabulatedVectorValueFvPatchVectorField&,
            const DimensionedField<vector, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchVectorField> clone
        (
            const DimensionedField<vector, volMesh>& iF
        ) const
        {
            return tmp<fvPatchVectorField>
            (
                new fixedTabulatedVectorValueFvPatchVectorField(*this, iF)
            );
        }


    // Member functions

        //- Return fileDir_
        std::string& fileDir()
        {
            return fileDir_;
        }
        //- Return nbPoint_
        unsigned int& nbPoint()
        {
            return nbPoint_;
        }
        //- Return wFactor_
        scalar& wFactor()
        {
            return wFactor_;
        }


        //- Update coefficients
        virtual void updateCoeffs();

        //- Write
        virtual void write(Ostream&) const;
};


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

} // End namespace Foam

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

#endif

// ************************************************************************* //
<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span><span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • fixedTabulatedVectorValueFvPatchScalarField.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
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

#include "fixedTabulatedVectorValueFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"


#include "nanoflann.hpp"
#include "kdtree.h"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


namespace Foam
{

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

fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(p, iF),
    fileDir_("."),
    nbPoint_(4),
    wFactor_(2)
{}


fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fixedTabulatedVectorValueFvPatchVectorField& ptf,
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
    fileDir_(ptf.fileDir_),
    nbPoint_(ptf.nbPoint_),
    wFactor_(ptf.wFactor_)
{}


fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    fileDir_(string(dict.lookup("fileDir"))),
    nbPoint_(dict.lookupOrDefault<scalar>("nbPoint",4)),
    wFactor_(dict.lookupOrDefault<scalar>("wFactor",2))
{

// check the existance of the file

    std::ifstream inFile(fileDir_,std::ios::in);  
 
    if(inFile) 
    {       
        inFile.close(); 
    }
    else 
    {
        FatalErrorIn("fixedTabulatedVectorValueFvPatchVectorField(dict)")
            << "path to the raw file is not correct"
            << abort(FatalError);
    }

    if( (nbPoint_ < 1) || (wFactor_ <= 0) ) 
    {
        FatalErrorIn("fixedTabulatedVectorValueFvPatchVectorField(dict)")
            << "the number of interpolation point is less than 1 or the weighting factor is not positive."
            << abort(FatalError);
    }

    evaluate();
}


fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fixedTabulatedVectorValueFvPatchVectorField& fcvpvf,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(fcvpvf, iF),
    fileDir_(fcvpvf.fileDir_),
    nbPoint_(fcvpvf.nbPoint_),
    wFactor_(fcvpvf.wFactor_)
{}


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

void fixedTabulatedVectorValueFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
// >> open the input data file
	std::ifstream inFile(fileDir_,std::ios::in);
        // count the lines in the file
	inFile.unsetf(std::ios_base::skipws);
    	unsigned nbLines = std::count(
        	std::istream_iterator<char>(inFile),
        	std::istream_iterator<char>(), 
        	'\n');
	nbLines=nbLines-2;
	//...skip two first lines (commentary)
	inFile.close();
	inFile.open(fileDir_);
	std::string line;  
        std::getline(inFile, line);  
        std::getline(inFile, line);  

	//Initialize a plain continous array for the data
        scalar vectorData[nbLines][6];
	unsigned int index_(0);

	PointCloud<scalar> cloud;
	//...for all line until the end, create a VectorField with the computed point
        while(std::getline(inFile, line))  // tant que l'on peut mettre la ligne dans "contenu"
        {
		std::istringstream buffer(line);
		std::string word_("");

	        buffer >> word_;
		scalar x_(atof(word_.c_str()));
	        buffer >> word_;
		scalar y_(atof(word_.c_str()));
	        buffer >> word_;
		scalar z_(atof(word_.c_str()));
	        buffer >> word_;
		scalar valx_(atof(word_.c_str()));
	        buffer >> word_;
		scalar valy_(atof(word_.c_str()));
	        buffer >> word_;
		scalar valz_(atof(word_.c_str()));

		updatePointCloud(cloud, x_, y_, z_); 

		vectorData[index_][0] = x_;
                vectorData[index_][1] = y_;
                vectorData[index_][2] = z_;
                vectorData[index_][3] = valx_;
		vectorData[index_][4] = valy_;
		vectorData[index_][5] = valz_;
		index_++;
	
        }

	
	// construct a kd-tree index:
	typedef nanoflann::KDTreeSingleIndexAdaptor<
		nanoflann::L2_Simple_Adaptor<scalar, PointCloud<scalar> > ,
		PointCloud<scalar>,
		3
		> my_kd_tree_t;

	my_kd_tree_t   index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(10) );
	index.buildIndex();

#if 0
	// Test resize of dataset and rebuild of index:
	cloud.pts.resize(cloud.pts.size()*0.5);
	index.buildIndex();
#endif


// >> get the face of the current patch & initialize a VectorField of the right size
	const vectorField& faceCenter = patch().Cf();
	vectorField values(faceCenter.size(), vector(0,0,0));

// >> compute the scalar value for all face centers
        forAll(faceCenter, f)
        {
// >>>>>>> get face center point
                const point faceCenterPoint = faceCenter[f];
// >>>>>>> find closest neighbourg
		// define the requested point
		const scalar query_pt[3] = { faceCenterPoint.x(), faceCenterPoint.y(), faceCenterPoint.z()};
		// find the 8 closest points
		size_t num_results = nbPoint_;
		std::vector<size_t>   ret_index(num_results);
		std::vector<scalar> out_dist_sqr(num_results);

		num_results = index.knnSearch(&query_pt[0], num_results, &ret_index[0], &out_dist_sqr[0]);
		
		// In case of less points in the tree than requested:
		ret_index.resize(num_results);
		out_dist_sqr.resize(num_results);

		// weight factors for PID
		std::vector<scalar> out_dist_weight(num_results);
		for (size_t i = 0; i < num_results; i++)
		{
			out_dist_weight[i]=1/std::pow(out_dist_sqr[i],wFactor_);
		}	
		// Distance inverse weighting methode (PID)
		scalar WxScalar(0);
		scalar WyScalar(0);
		scalar WzScalar(0);
		scalar Wtot(0);
		for (size_t i = 0; i < num_results; i++)
		{
			WxScalar+=out_dist_weight[i]*vectorData[ret_index[i]][3];
			WyScalar+=out_dist_weight[i]*vectorData[ret_index[i]][4];
			WzScalar+=out_dist_weight[i]*vectorData[ret_index[i]][5];
			Wtot+=out_dist_weight[i];
		}
		values[f]=vector(WxScalar/Wtot, WyScalar/Wtot, WzScalar/Wtot); ;


        }



    // Calculate local 1-D coordinate for the parabolic profile
   // VectorField coord = 0; //2*((c - ctr) & y_)/((bb.max() - bb.min()) & y_);

   vectorField::operator=(values);
}


// Write
void fixedTabulatedVectorValueFvPatchVectorField::write(Ostream& os) const
{
    fvPatchVectorField::write(os);
    os.writeKeyword("fileDir")
        << fileDir_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


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

makePatchTypeField(fvPatchVectorField, fixedTabulatedVectorValueFvPatchVectorField);

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

} // End namespace Foam

// ************************************************************************* //
<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span><span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

 

Non classé

OpenFOAM – fixedTabulatedScalarValue

Alexis Sauvageon

Cette condition limite permet d’imposer au bord du domaine fluide des valeurs scalaires stockées dans une table. Les valeurs scalaires doivent être associées à un point définit par ses coordonnées.

xyzepsilon
-12-4200.00800494
-12-420.2740780.00730994
-12-41.60640.2740780.00792009
⋮⋮⋮⋮

Le programme utilise un algorithme des plus proches voisins pour associer à chaque centre de faces du patch sur lequel est imposée la condition limite une valeur d’entrée calculée par la méthode de la pondération inverse à la distance.

Les développements font appel à la librarie nanoflann accessible en cliquant sur le lien suivant. Cette librarie ne requiert aucune compilation préalable, toutefois les liens vers les headers doivent être correctement définis pour pouvoir compiler les sources de fixedTabularScalarValue.

Dans la suite, nous fournissons les sources du .h et du .cpp ainsi que des deux fichiers options et files nécessaires à la compilation. Les chemins sont susceptibles de changer en fonction de l’environnement de travail.

 

  • options
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(FOAM_DEV)/libraries/nanoflann/include \
-I$(FOAM_DEV)/libraries/utils/include
EXE_LIBS =<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • files
fvPatchFields = ./
derivedFvPatchFields = $(fvPatchFields)/derived

$(derivedFvPatchFields)/fixedTabulatedScalarValue/fixedTabulatedScalarValueFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/libArepCustomPatch

 

  • fixedTabulatedScalarValueFvPatchScalarField.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Class
    fixedTabulatedScalarValueFvPatchScalarField

Description
    Compute boundaries conditions from a field stored in table (format .raw)
	[X Y Z Scalar]
    the two first lines of the table are not interpreted.
    The interpollated value is computed thanks to a PID method.

    \heading Patch usage

    \table
        Property    | Description             	 	| Required   	| Default value
        fileDir     | path to the input file  		| yes        	| 
        nbPoint     | number of point for PID 	        | no       	| 4
        wFactor     | p factor for PID      		| no        	| 2.0

    \endtable

    Example of the boundary condition specification:
    \verbatim
    myPatch
    {
        type            fixedTabulatedScalarValue;
        fileDir         "constant/inputTableData/p_boxValue.raw";
        nbPoint              8;
        wFactor            2.0;

    }
    \endverbatim


SourceFiles
    fixedTabulatedScalarValueFvPatchScalarField.C

Author
    Alexis Sauvageon.  All rights reserved

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

#ifndef fixedTabulatedScalarValueFvPatchScalarField_H
#define fixedTabulatedScalarValueFvPatchScalarField_H

#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"


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

namespace Foam
{

/*---------------------------------------------------------------------------*\
              Class fixedTabulatedScalarValueFvPatchField Declaration
\*---------------------------------------------------------------------------*/

class fixedTabulatedScalarValueFvPatchScalarField
:
    public fixedValueFvPatchScalarField
{
    // Private data

        //- Peak velocity magnitude
        std::string fileDir_;
	unsigned int nbPoint_;
	scalar wFactor_;

public:

    //- Runtime type information
    TypeName("fixedTabulatedScalarValue");


    // Constructors

        //- Construct from patch and internal field
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const dictionary&
        );

        //- Construct by mapping given fixedTabulatedScalarValueFvPatchScalarField
        //  onto a new patch
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fixedTabulatedScalarValueFvPatchScalarField&,
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchScalarField> clone() const
        {
            return tmp<fvPatchScalarField>
            (
                new fixedTabulatedScalarValueFvPatchScalarField(*this)
            );
        }

        //- Construct as copy setting internal field reference
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fixedTabulatedScalarValueFvPatchScalarField&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchScalarField> clone
        (
            const DimensionedField<scalar, volMesh>& iF
        ) const
        {
            return tmp<fvPatchScalarField>
            (
                new fixedTabulatedScalarValueFvPatchScalarField(*this, iF)
            );
        }


    // Member functions

        //- Return fileDir_
        std::string& fileDir()
        {
            return fileDir_;
        }
        //- Return nbPoint_
        unsigned int& nbPoint()
        {
            return nbPoint_;
        }
        //- Return wFactor_
        scalar& wFactor()
        {
            return wFactor_;
        }


        //- Update coefficients
        virtual void updateCoeffs();

        //- Write
        virtual void write(Ostream&) const;
};


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

} // End namespace Foam

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

#endif

// ************************************************************************* //
<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • fixedTabulatedScalarValueFvPatchScalarField.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
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

#include "fixedTabulatedScalarValueFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include <fstream>
#include <sstream>
#include <iterator>

#include "nanoflann.hpp"
#include "kdtree.h"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


namespace Foam
{

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

fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedValueFvPatchScalarField(p, iF),
    fileDir_("."),
    nbPoint_(4),
    wFactor_(2)
{}


fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fixedTabulatedScalarValueFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
    fileDir_(ptf.fileDir_),
    nbPoint_(ptf.nbPoint_),
    wFactor_(ptf.wFactor_)
{}


fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchScalarField(p, iF),
    fileDir_(string(dict.lookup("fileDir"))),
    nbPoint_(dict.lookupOrDefault<scalar>("nbPoint",4)),
    wFactor_(dict.lookupOrDefault<scalar>("wFactor",2))
{

// check the existance of the file

    std::ifstream inFile(fileDir_,std::ios::in);  
 
    if(inFile) 
    {       
        inFile.close(); 
    }
    else 
    {
        FatalErrorIn("fixedTabulatedScalarValueFvPatchScalarField(dict)")
            << "path to the raw file is not correct"
            << abort(FatalError);
    }

    if( (nbPoint_ < 1) || (wFactor_ <= 0) ) 
    {
        FatalErrorIn("fixedTabulatedScalarValueFvPatchScalarField(dict)")
            << "the number of interpolation point is less than 1 or the weighting factor is not positive."
            << abort(FatalError);
    }

    evaluate();
}


fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fixedTabulatedScalarValueFvPatchScalarField& fcvpvf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedValueFvPatchScalarField(fcvpvf, iF),
    fileDir_(fcvpvf.fileDir_),
    nbPoint_(fcvpvf.nbPoint_),
    wFactor_(fcvpvf.wFactor_)
{}


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

void fixedTabulatedScalarValueFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
// >> open the input data file
	std::ifstream inFile(fileDir_,std::ios::in);
        // count the lines in the file
	inFile.unsetf(std::ios_base::skipws);
    	unsigned nbLines = std::count(
        	std::istream_iterator<char>(inFile),
        	std::istream_iterator<char>(), 
        	'\n');
	nbLines=nbLines-2;
	//...skip two first lines (commentary)
	inFile.close();
	inFile.open(fileDir_);
	std::string line;  
        std::getline(inFile, line);  
        std::getline(inFile, line);  

	//Initialize a plain continous array for the data
        scalar vectorData[nbLines][4];
	unsigned int index_(0);

	PointCloud<scalar> cloud;
	//...for all line until the end, create a scalarField with the computed point
        while(std::getline(inFile, line))  // tant que l'on peut mettre la ligne dans "contenu"
        {
		std::istringstream buffer(line);
		std::string word_("");

	        buffer >> word_;
		scalar x_(atof(word_.c_str()));
	        buffer >> word_;
		scalar y_(atof(word_.c_str()));
	        buffer >> word_;
		scalar z_(atof(word_.c_str()));
	        buffer >> word_;
		scalar val_(atof(word_.c_str()));

		updatePointCloud(cloud, x_, y_, z_); 

		vectorData[index_][0] = x_;
                vectorData[index_][1] = y_;
                vectorData[index_][2] = z_;
                vectorData[index_][3] = val_;
		index_++;
	
        }

	
	// construct a kd-tree index:
	typedef nanoflann::KDTreeSingleIndexAdaptor<
		nanoflann::L2_Simple_Adaptor<scalar, PointCloud<scalar> > ,
		PointCloud<scalar>,
		3
		> my_kd_tree_t;

	my_kd_tree_t   index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(10) );
	index.buildIndex();

#if 0
	// Test resize of dataset and rebuild of index:
	cloud.pts.resize(cloud.pts.size()*0.5);
	index.buildIndex();
#endif


// >> get the face of the current patch & initialize a scalarField of the right size
	const vectorField& faceCenter = patch().Cf();
	scalarField values(faceCenter.size(), scalar(0));

// >> compute the scalar value for all face centers
        forAll(faceCenter, f)
        {
// >>>>>>> get face center point
                const point faceCenterPoint = faceCenter[f];
// >>>>>>> find closest neighbourg
		// define the requested point
		const scalar query_pt[3] = { faceCenterPoint.x(), faceCenterPoint.y(), faceCenterPoint.z()};
		// find the 8 closest points
		size_t num_results = nbPoint_;
		std::vector<size_t>   ret_index(num_results);
		std::vector<scalar> out_dist_sqr(num_results);

		num_results = index.knnSearch(&query_pt[0], num_results, &ret_index[0], &out_dist_sqr[0]);
		
		// In case of less points in the tree than requested:
		ret_index.resize(num_results);
		out_dist_sqr.resize(num_results);

		// weight factors for PID
		std::vector<scalar> out_dist_weight(num_results);
		for (size_t i = 0; i < num_results; i++)
		{
			out_dist_weight[i]=1/std::pow(out_dist_sqr[i],wFactor_);
		}	
		// Distance inverse weighting methode (PID)
		scalar WxScalar(0);
		scalar Wtot(0);
		for (size_t i = 0; i < num_results; i++)
		{
			WxScalar+=out_dist_weight[i]*vectorData[ret_index[i]][3];
			Wtot+=out_dist_weight[i];
		}
		values[f]=WxScalar/Wtot;

		// print the index of the closest point and their distance
		/*cout << faceCenterPoint.x() << "    " << faceCenterPoint.y() << "    " << faceCenterPoint.z()  << std::endl;
		cout << values[f] << std::endl;
		cout << std::endl;

		for (size_t i = 0; i < num_results; i++)
		{
		cout << vectorData[ret_index[i]][0]  << "    " << vectorData[ret_index[i]][1] << "    " << vectorData[ret_index[i]][2] << std::endl;
		cout << "dist : " << out_dist_sqr[i] << " value : " << vectorData[ret_index[i]][4]  << std::endl;
		cout << "weight : " << out_dist_weight[i]   << std::endl;
		}*/

        }



    // Calculate local 1-D coordinate for the parabolic profile
   // scalarField coord = 0; //2*((c - ctr) & y_)/((bb.max() - bb.min()) & y_);

   scalarField::operator=(values);
}


// Write
void fixedTabulatedScalarValueFvPatchScalarField::write(Ostream& os) const
{
    fvPatchScalarField::write(os);
    os.writeKeyword("fileDir")
        << fileDir_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


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

makePatchTypeField(fvPatchScalarField, fixedTabulatedScalarValueFvPatchScalarField);

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

} // End namespace Foam

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