errorEstimate

In this article, we propose a slightly modified version of an error estimator proposed by the foam-extend fork.

This is a residual error estimator used to estimate the discretisation error due to the finite volume method on the elements of a 2D or 3D mesh. For more information on error estimators, see the work of [H. Jsak].( https://spiral.imperial.ac.uk/bitstream/10044/1/8335/1/Hrvoje_Jasak-1996-PhD-Thesis.pdf ) or in the field of continuum mechanics at Code_Aster documentation .

This version of the error estimator, recoded to be compatible with OpenFoam v5.0 and available here, is indicative and has yet to be validated.

• files

LIB = $(FOAM_LIBBIN)/liberrorEstimation

• option

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

LIB_LIBS = \
    -lfiniteVolume \
    -lmeshTools

• errorEstimate.H

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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

  Show code

  /*---------------------------------------------------------------------------*\
    =========                 |
    \\      /  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
  
  // ************************************************************************* //
  
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //