template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> compareLR
(
  const GeometricField<Type, fvsPatchField, surfaceMesh>& sfL,
  const GeometricField<Type, fvsPatchField, surfaceMesh>& sfR,
  const surfaceScalarField& select,
  const scalar thres
)
{
  const fvMesh& mesh = sfL.mesh();
  tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsf
  (
    new GeometricField<Type, fvsPatchField, surfaceMesh>
        (
            IOobject
            (
                "select("+sfL.name()+","+sfR.name()+')',
                sfL.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            sfL.dimensions()
         )
    );
  GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();

  forAll(select, facei)
    {
       sf[facei] = (select[facei] >= thres ? sfL[facei] : sfR[facei]);
    }
  
  forAll(select.boundaryField(), patchi)
    {
        Field<Type>& psf = sf.boundaryFieldRef()[patchi];
        const Field<Type>& psfL = sfL.boundaryField()[patchi];
        const Field<Type>& psfR = sfR.boundaryField()[patchi];
        const scalarField& pSel = select.boundaryField()[patchi];
        forAll(pSel, facei)
        {
            psf[facei] = (pSel[facei] >= thres ? psfL[facei] : psfR[facei]);   
        }
    }

  return tsf;
}

template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > M1
(
  const GeometricField<Type, fvsPatchField, surfaceMesh>& Ma,
  const GeometricField<Type, fvsPatchField, surfaceMesh>& sgn
)
{
  const fvMesh& mesh = Ma.mesh();
  tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf
  (
    new GeometricField<Type, fvsPatchField, surfaceMesh>
        (
            IOobject
            (
                "Ma1"+sgn.name(),
                Ma.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            Ma.dimensions()
         )
    );
  GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();

  forAll(Ma, facei)
     {
       sf[facei] = 0.5*(Ma[facei] + sgn[facei]*mag(Ma[facei]));
     }
 
  forAll(Ma.boundaryField(), patchi)
     {
       Field<Type>& psf = sf.boundaryFieldRef()[patchi];
       const Field<Type>& pMa = Ma.boundaryField()[patchi];
       const Field<Type>& psgn = sgn.boundaryField()[patchi];
       forAll(psf, facei)
          {
            psf[facei] = 0.5*(pMa[facei] + psgn[facei]*mag(pMa[facei]));
           }
      }

   return tsf;
}

template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > M2
(
  const GeometricField<Type, fvsPatchField, surfaceMesh>& Ma,
  const GeometricField<Type, fvsPatchField, surfaceMesh>& sgn
)
{
  const fvMesh& mesh = Ma.mesh();
  tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf
  (
    new GeometricField<Type, fvsPatchField, surfaceMesh>
        (
            IOobject
            (
                "Ma2"+sgn.name(),
                Ma.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            Ma.dimensions()
         )
    );
  GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();

  forAll(Ma, facei)
     {
       sf[facei] = sgn[facei]*0.25*(Ma[facei] + sgn[facei])*(Ma[facei] + sgn[facei]);
     }
 
  forAll(Ma.boundaryField(), patchi)
     {
       Field<Type>& psf = sf.boundaryFieldRef()[patchi];
       const Field<Type>& pMa = Ma.boundaryField()[patchi];
       const Field<Type>& psgn = sgn.boundaryField()[patchi];
       forAll(psf, facei)
          {
            psf[facei] = psgn[facei]*0.25*(pMa[facei] + psgn[facei])*(pMa[facei] + psgn[facei]);
          }
      }

   return tsf;
}

template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > M4
(
  const GeometricField<Type, fvsPatchField, surfaceMesh>& Ma,
  const GeometricField<Type, fvsPatchField, surfaceMesh>& sgn,
  scalar beta
)
{
  const fvMesh& mesh = Ma.mesh();
  tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf
  (
    new GeometricField<Type, fvsPatchField, surfaceMesh>
        (
            IOobject
            (
                "Ma4"+sgn.name(),
                Ma.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            Ma.dimensions()
         )
    );
  const GeometricField<Type, fvsPatchField, surfaceMesh>& rsgn = -sgn;
  GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();
  GeometricField<Type, fvsPatchField, surfaceMesh> Ma2pos = M2(Ma, sgn);
  GeometricField<Type, fvsPatchField, surfaceMesh> Ma2neg = M2(Ma, rsgn);
  GeometricField<Type, fvsPatchField, surfaceMesh> Ma1 = M1(Ma, sgn);
  forAll(Ma, facei)
     {
       if (mag(Ma[facei]) < 1.)
        {
         sf[facei] = Ma2pos[facei]*(1.-sgn[facei]*16.*beta*Ma2neg[facei]);
        }
       else
        {
         sf[facei] = Ma1[facei];
        }
     }
    
  forAll(Ma.boundaryField(), patchi)
     {
       Field<Type>& psf = sf.boundaryFieldRef()[patchi];
       const Field<Type>& pMa = Ma.boundaryField()[patchi];
       const Field<Type>& pMa2pos = Ma2pos.boundaryField()[patchi];
       const Field<Type>& pMa2neg = Ma2neg.boundaryField()[patchi];
       const Field<Type>& pMa1 = Ma1.boundaryField()[patchi];
       const Field<Type>& psgn = sgn.boundaryField()[patchi];
       forAll(psf, facei)
          {
            if (mag(pMa[facei])<1)
             {
              psf[facei] = pMa2pos[facei]*(1.-psgn[facei]*16.*beta*pMa2neg[facei]);
             }
            else
             {
              psf[facei] = pMa1[facei];
             }
          }
      }

   return tsf;
}

template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > p5
(
  const GeometricField<Type, fvsPatchField, surfaceMesh>& Ma,
  const GeometricField<Type, fvsPatchField, surfaceMesh>& sgn,
  const GeometricField<Type, fvsPatchField, surfaceMesh>& alpha
)
{
  const fvMesh& mesh = Ma.mesh();
  tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf
  (
    new GeometricField<Type, fvsPatchField, surfaceMesh>
        (
            IOobject
            (
                "p5"+sgn.name(),
                Ma.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            Ma.dimensions()
         )
    );
  const GeometricField<Type, fvsPatchField, surfaceMesh>& rsgn = -sgn;
  GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();
  GeometricField<Type, fvsPatchField, surfaceMesh> Ma2pos = M2(Ma, sgn);
  GeometricField<Type, fvsPatchField, surfaceMesh> Ma2neg = M2(Ma, rsgn);
  GeometricField<Type, fvsPatchField, surfaceMesh> Ma1 = M1(Ma, sgn);

  forAll(Ma, facei)
     {
       if (mag(Ma[facei]) < 1.)
        {
         sf[facei] = Ma2pos[facei]*( (sgn[facei]*2. - Ma[facei]) - sgn[facei]*16.*alpha[facei]*Ma[facei]*Ma2neg[facei]);
        }
       else
        {
         sf[facei] = 1./Ma[facei]*Ma1[facei];
        }
     }
  
  forAll(Ma.boundaryField(), patchi)
     {
       Field<Type>& psf = sf.boundaryFieldRef()[patchi];
       const Field<Type>& pMa = Ma.boundaryField()[patchi];
       const Field<Type>& pMa2pos = Ma2pos.boundaryField()[patchi];
       const Field<Type>& pMa2neg = Ma2neg.boundaryField()[patchi];
       const Field<Type>& pMa1 = Ma1.boundaryField()[patchi];
       const Field<Type>& psgn = sgn.boundaryField()[patchi];
       const Field<Type>& palpha = alpha.boundaryField()[patchi];
       forAll(psf, facei)
          {
            if (mag(pMa[facei]) < 1.)
             {
              psf[facei] = pMa2pos[facei]*((psgn[facei]*2. - pMa[facei])-psgn[facei]*16.*palpha[facei]*pMa[facei]*pMa2neg[facei]);
             }
            else
             {
              psf[facei] = 1./pMa[facei]*pMa1[facei];
             }
          }
      }

   return tsf;
}




















