Convection.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  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, 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: Convection.hpp,v 1.17 2007/04/22 21:46:10 delpinux Exp $


#ifndef CONVECTION_HPP
#define CONVECTION_HPP

#include <ScalarFunctionBase.hpp>
#include <FieldOfScalarFunction.hpp>

#include <Structured3DMesh.hpp>
#include <MeshOfTetrahedra.hpp>

#include <ConnectivityBuilder.hpp>

#include <ConformTransformation.hpp>

#include <limits>
template <typename MeshType>
struct Convection
  : public ScalarFunctionBase
{
};

template <>
class Convection<Structured3DMesh>
  : public ScalarFunctionBase
{
private:
  std::ostream& __put(std::ostream& os) const
  {
    os << "convect(" << __u << ","
       << __deltaT << ',' << __phi << ')';
    return os;
  }

  const real_t __deltaT;

  const Structured3DMesh& __mesh;

  const FieldOfScalarFunction& __u;

  const ScalarFunctionBase& __phi;

  Convection(const Convection<Structured3DMesh>& C);

public:
  real_t operator()(const TinyVector<3>& X) const
  {
    Structured3DMesh::const_iterator h = __mesh.find(X[0], X[1], X[2]);
    if (h.end()) {
      return 0; // we are outside the mesh
    }

    const ScalarFunctionBase& u0 = *__u.function(0);
    const ScalarFunctionBase& u1 = *__u.function(1);
    const ScalarFunctionBase& u2 = *__u.function(2);

    const Connectivity<Structured3DMesh>& connectivity = __mesh.connectivity();
    if (not(connectivity.hasCellToCells())) {
      ConnectivityBuilder<Structured3DMesh> c(__mesh);
      c.generates(Connectivity<Structured3DMesh>::CellToCells);
    }

    TinyVector<3, real_t> X0 = X;

    TinyVector<3, real_t> Xhat0;
    {
      ConformTransformationQ1CartesianHexahedron T0(*h);
      T0.invertT(X0, Xhat0);
    }

    real_t dt0 = __deltaT;

    const size_t faceNumberConverter[6] = {4,2,1,3,0,5};

    do {
      const CartesianHexahedron& H = (*h);
      ConformTransformationQ1CartesianHexahedron T(H);

      T.value(Xhat0,X0);

      TinyVector<3, real_t> X1;
      X1[0] = X0[0]-u0(X0)*dt0;
      X1[1] = X0[1]-u1(X0)*dt0;
      X1[2] = X0[2]-u2(X0)*dt0;

      TinyVector<3, real_t> Xhat1;
      T.invertT(X1, Xhat1);

      if((   Xhat1[0] >= 0)
         and(Xhat1[1] >= 0)
         and(Xhat1[2] >= 0)
         and(Xhat1[0] <= 1)
         and(Xhat1[1] <= 1)
         and(Xhat1[2] <= 1)) {
        dt0 = 0;
        X0 = X1;
      } else {

        const TinyVector<3, real_t> deltaX = Xhat1-Xhat0;
        
        // Get the planes which are required to compute intersection
        TinyVector<3, real_t> x0 = 0;
        for (size_t i=0; i<3; ++i) {
          if (deltaX[i]>0) x0[i] = 1;
        }

        // Computes the linear abcisse required to reach faces
        TinyVector<3, real_t> coef = std::numeric_limits<real_t>::max();
        for (size_t i=0; i<3; ++i) {
          if (std::abs(deltaX[i])>0) {
            coef[i] = std::abs((x0[i]-Xhat0[i])/deltaX[i]);
          }
        }

        size_t minimumCoefficentNumber = 0;
        for (size_t i=1; i<3; ++i) {
          minimumCoefficentNumber
            = (coef[minimumCoefficentNumber]<coef[i])?minimumCoefficentNumber:i;
        }

        const size_t outGoingFace
          = faceNumberConverter[2*minimumCoefficentNumber + ((deltaX[minimumCoefficentNumber]>0)?1:0)];

        const real_t dt = dt0*coef[minimumCoefficentNumber];

        Xhat0 += coef[minimumCoefficentNumber] * deltaX;
        dt0-=dt;

        h = connectivity.cells(*h)[outGoingFace];
        if (h.end()) {
          T.value(Xhat0,X0);
          dt0=0;
        } // takes the value on the border when going outside

        // X0 is not updated. We compute the value of Xhat0 in the next cell
        // This make the periodic boundary conditions valid
        Xhat0[minimumCoefficentNumber] = 1-x0[minimumCoefficentNumber];
      }
    } while (dt0>0);
    return __phi(X0);
  }

  bool canBeSimplified() const
  {
    return false;
  }

  Convection(const ScalarFunctionBase& phi,
             const FieldOfScalarFunction& u,
             const real_t& dt,
             const Structured3DMesh& M)
    : ScalarFunctionBase(ScalarFunctionBase::convection),
      __deltaT(dt),
      __mesh(M),
      __u(u),
      __phi(phi)
  {
    if (__u.numberOfComponents() != 3) {
      throw ErrorHandler(__FILE__,__LINE__,
                         "Convection needs a 3 component field, you provided "+stringify(__u),
                         ErrorHandler::normal);
    }
  }

  ~Convection()
  {
    ;
  }
};

template <>
class Convection<MeshOfTetrahedra>
  : public ScalarFunctionBase
{
private:
  std::ostream& __put(std::ostream& os) const
  {
    os << "convect(" << __u << ","
       << __deltaT << ',' << __phi << ')';
    return os;
  }

  const real_t __deltaT;

  const MeshOfTetrahedra& __mesh;

  const FieldOfScalarFunction& __u;

  const ScalarFunctionBase& __phi;
public:
  real_t operator()(const TinyVector<3>& X) const
  {
    const real_t epsilon = 1E-6;
    MeshOfTetrahedra::const_iterator t = __mesh.find(X[0], X[1], X[2]);
    if (t.end()) return 0; // we are outside the mesh

    const ScalarFunctionBase& u0 = *__u.function(0);
    const ScalarFunctionBase& u1 = *__u.function(1);
    const ScalarFunctionBase& u2 = *__u.function(2);

    const Connectivity<MeshOfTetrahedra>& connectivity = __mesh.connectivity();
    if (not(connectivity.hasCellToCells())) {
      ConnectivityBuilder<MeshOfTetrahedra> c(__mesh);
      c.generates(Connectivity<MeshOfTetrahedra>::CellToCells);
    }

    TinyVector<3, real_t> X0 = X;

    real_t dt0 = __deltaT;
    do {
      TinyVector<3, real_t> X1;
      X1[0] = X0[0]-u0(X0)*dt0;
      X1[1] = X0[1]-u1(X0)*dt0;
      X1[2] = X0[2]-u2(X0)*dt0;

      TinyVector<4, real_t> lambda1;
      const Tetrahedron& T = (*t);
      T.getBarycentricCoordinates(X1, lambda1);
      if((lambda1[0]>-epsilon)
         and(lambda1[1]>-epsilon)
         and(lambda1[2]>-epsilon)
         and(lambda1[3]>-epsilon)) {
        dt0 = 0;
        X0 = X1;
      } else {
        TinyVector<4, real_t> lambda0;
        T.getBarycentricCoordinates(X0, lambda0);

        TinyVector<4, real_t> deltaLambda = lambda1-lambda0;
        real_t coeff = std::numeric_limits<real_t>::max();
        size_t outGoingFace=std::numeric_limits<size_t>::max();
        for (size_t i=0; i<4; ++i) {
          if (lambda1[i] <= -epsilon) {
            real_t tmp = -lambda0[i]/deltaLambda[i];
            if ((std::abs(tmp) < std::abs(coeff))) {
              coeff = tmp;
              outGoingFace = i;
            }
          }
        }
        if (not(std::abs(lambda0[outGoingFace]) < epsilon)) { // we are not going back to the previous element
          const real_t dt = dt0*coeff;
          dt0-=dt;
          X0 += coeff * (X1-X0);
          T.getBarycentricCoordinates(X0, lambda0);
        }

        t = connectivity.cells(*t)[outGoingFace];
        if (t.end()) { dt0=0; } // takes the value on the border when going outside
      }
    } while (dt0>0);
    return __phi(X0);
  }

  bool canBeSimplified() const
  {
    return false;
  }

  Convection(const Convection<MeshOfTetrahedra>& C)
    : ScalarFunctionBase(ScalarFunctionBase::convection),
      __deltaT(C.__deltaT),
      __mesh(C.__mesh),
      __u(C.__u),
      __phi(C.__phi)
  {
    ;
  }

  Convection(const ScalarFunctionBase& phi,
             const FieldOfScalarFunction& u,
             const real_t& dt,
             const MeshOfTetrahedra& M)
    : ScalarFunctionBase(ScalarFunctionBase::convection),
      __deltaT(dt),
      __mesh(M),
      __u(u),
      __phi(phi)
  {
    ;
  }

  ~Convection()
  {
    ;
  }
};
template <>
class Convection<MeshOfHexahedra>
  : public ScalarFunctionBase
{
private:
  std::ostream& __put(std::ostream& os) const
  {
    os << "convect(" << __u << ","
       << __deltaT << ',' << __phi << ')';
    return os;
  }

  const real_t __deltaT;

  const MeshOfHexahedra& __mesh;

  const FieldOfScalarFunction& __u;

  const ScalarFunctionBase& __phi;

public:
  real_t operator()(const TinyVector<3>& X) const
  {
    MeshOfHexahedra::const_iterator h = __mesh.find(X[0], X[1], X[2]);
    if (h.end()) return 0; // we are outside the mesh

    const ScalarFunctionBase& u0 = *__u.function(0);
    const ScalarFunctionBase& u1 = *__u.function(1);
    const ScalarFunctionBase& u2 = *__u.function(2);

    const Connectivity<MeshOfHexahedra>& connectivity = __mesh.connectivity();
    if (not(connectivity.hasCellToCells())) {
      ConnectivityBuilder<MeshOfHexahedra> c(__mesh);
      c.generates(Connectivity<MeshOfHexahedra>::CellToCells);
    }

    TinyVector<3, real_t> X0 = X;

    TinyVector<3, real_t> Xhat0;
    {
      ConformTransformationQ1Hexahedron T0(*h);
      T0.invertT(X0, Xhat0);
    }

    real_t dt0 = __deltaT;

    const size_t faceNumberConverter[6] = {4,2,1,3,0,5};

    do {
      const Hexahedron& H = (*h);
      ConformTransformationQ1Hexahedron T(H);

      T.value(Xhat0,X0);

      TinyVector<3, real_t> X1;
      X1[0] = X0[0]-u0(X0)*dt0;
      X1[1] = X0[1]-u1(X0)*dt0;
      X1[2] = X0[2]-u2(X0)*dt0;

      TinyVector<3, real_t> Xhat1;
      T.invertT(X1, Xhat1);

      if((   Xhat1[0] >= 0)
         and(Xhat1[1] >= 0)
         and(Xhat1[2] >= 0)
         and(Xhat1[0] <= 1)
         and(Xhat1[1] <= 1)
         and(Xhat1[2] <= 1)) {
        dt0 = 0;
        X0 = X1;
      } else {

        const TinyVector<3, real_t> deltaX = Xhat1-Xhat0;
        
        // Get the planes which are required to compute intersection
        TinyVector<3, real_t> x0 = 0;
        for (size_t i=0; i<3; ++i) {
          if (deltaX[i]>0) x0[i] = 1;
        }

        // Computes the linear abcisse required to reach faces
        TinyVector<3, real_t> coef = std::numeric_limits<real_t>::max();
        for (size_t i=0; i<3; ++i) {
          if (std::abs(deltaX[i])>0) {
            coef[i] = std::abs((x0[i]-Xhat0[i])/deltaX[i]);
          }
        }

        size_t minimumCoefficentNumber = 0;
        for (size_t i=1; i<3; ++i) {
          minimumCoefficentNumber
            = (coef[minimumCoefficentNumber]<coef[i])?minimumCoefficentNumber:i;
        }

        const size_t outGoingFace
          = faceNumberConverter[2*minimumCoefficentNumber + ((deltaX[minimumCoefficentNumber]>0)?1:0)];

        const real_t dt = dt0*coef[minimumCoefficentNumber];

        Xhat0 += coef[minimumCoefficentNumber] * deltaX;
        dt0-=dt;

        h = connectivity.cells(*h)[outGoingFace];
        if (h.end()) {
          T.value(Xhat0,X0);
          dt0=0;
        } // takes the value on the border when going outside

        // X0 is not updated. We compute the value of Xhat0 in the next cell
        // This make the periodic boundary conditions valid
        Xhat0[minimumCoefficentNumber] = 1-x0[minimumCoefficentNumber];
      }
    } while (dt0>0);

    return __phi(X0);
  }

  bool canBeSimplified() const
  {
    return false;
  }

  Convection(const Convection<MeshOfHexahedra>& C)
    : ScalarFunctionBase(ScalarFunctionBase::convection),
      __deltaT(C.__deltaT),
      __mesh(C.__mesh),
      __u(C.__u),
      __phi(C.__phi)
  {
    ;
  }

  Convection(const ScalarFunctionBase& phi,
             const FieldOfScalarFunction& u,
             const real_t& dt,
             const MeshOfHexahedra& M)
    : ScalarFunctionBase(ScalarFunctionBase::convection),
      __deltaT(dt),
      __mesh(M),
      __u(u),
      __phi(phi)
  {
    ;
  }

  ~Convection()
  {
    ;
  }
};

#endif // CONVECTION_HPP

---