text: #ifndef COOLFluiD_Numerics_FiniteVolume_LeastSquareP1PolyRec3D_hh
  #define COOLFluiD_Numerics_FiniteVolume_LeastSquareP1PolyRec3D_hh
  
  //////////////////////////////////////////////////////////////////////////////
  
  #include "FiniteVolume/FVMCC_PolyRec.hh"
  
  #ifdef CF_HAVE_CUDA
  #include "FiniteVolume/FluxData.hh"
  #include "FiniteVolume/KernelData.hh"
  #include "FiniteVolume/CellData.hh"
  #include "Common/CUDA/CudaEnv.hh"
  #include "Framework/SubSystemStatus.hh"
  #endif
  
  //////////////////////////////////////////////////////////////////////////////
  
  namespace COOLFluiD {
  
    namespace Numerics {
  
      namespace FiniteVolume {
  
  //////////////////////////////////////////////////////////////////////////////
  
  /**
   * This class implements a least square polynomial reconstructor in 3D for FVM
   *
   * @author Mehmet Sarp Yalim
   * @author Andrea Lani
   */
  class LeastSquareP1PolyRec3D : public FVMCC_PolyRec {
  public:
  
  #ifdef CF_HAVE_CUDA
    /// nested class defining local options
    template <typename P = NOTYPE>
    class DeviceConfigOptions {
    public:
      /// constructor
      HOST_DEVICE DeviceConfigOptions() {}
      
      /// destructor
      HOST_DEVICE ~DeviceConfigOptions() {}
      
      /// initialize with another object of the same kind
      HOST_DEVICE void init(DeviceConfigOptions<P> *const in) 
      {
        limitIter      = in->limitIter;
        currIter       = in->currIter; 
        limitRes       = in->limitRes; 
        currRes        = in->currRes;  
        gradientFactor = in->gradientFactor;
        freezeLimiter  = in->freezeLimiter;
      }
      
      CFuint limitIter;      /// iteration at which limiter is activated
      CFuint currIter;       /// current iteration
      CFreal limitRes;       /// residual at which limiter is activated
      CFreal currRes;        /// current residual
      CFreal gradientFactor; /// gradient factor within [0,1]
      CFuint freezeLimiter;    /// flag controlling the freezing of the limiter
    };
    
    /// nested class defining a functor
    template <typename PHYS>
    class DeviceFunc {
    public:
      typedef LeastSquareP1PolyRec3D BASE;
  
      /// constructor taking the options as argument
      HOST_DEVICE DeviceFunc(DeviceConfigOptions<NOTYPE>* dco) : m_dco(dco) {}
      
      /// compute the cell-based gradients
      HOST_DEVICE void computeGradients(const CFreal* state, const CFreal* node,
  				      KernelData<CFreal> *const kd,
  				      CellData::Itr* cell); 
      
      /// extrapolate the solution on the given face
      HOST_DEVICE void extrapolateOnFace(FluxData<PHYS>* currFd, const CFreal* coord, const CFreal* uX, 
  				       const CFreal* uY, const CFreal* uZ, const CFreal* limiter);
      
      /// extrapolate the solution on the given face for a given variable
      HOST_DEVICE inline void extrapolateOnFace(const CFuint iVar,
  					      FluxData<PHYS>* currFd, const CFreal* coord, const CFreal* uX, 
  					      const CFreal* uY, const CFreal* uZ, const CFreal* limiter);
      
      /// compute the distance between two points
      template <CFuint DIM, typename T1, typename T2>
      HOST_DEVICE inline CFreal getDistance(const T1* n1, const T2* n2)
      {
        CFreal dist = 0.;
        for (size_t i = 0; i < DIM; ++i) {
  	const CFreal diff = n1[i] - n2[i]; dist += diff*diff;
        }
        return sqrt(dist);
      }
      
    private:
      DeviceConfigOptions<NOTYPE>* m_dco;
    };
    
    /// copy the local configuration options to the device
    void copyConfigOptionsToDevice(DeviceConfigOptions<NOTYPE>* dco) 
    {
      // consider to copy to constant memory
      CudaEnv::copyHost2Dev(&dco->limitIter, &_limitIter, 1);
      CudaEnv::copyHost2Dev(&dco->limitRes, &_limitRes, 1);
      CudaEnv::copyHost2Dev(&dco->gradientFactor, &_gradientFactor, 1);
      
      CFuint flag = (_freezeLimiter) ? 1 : 0;
      CudaEnv::copyHost2Dev(&dco->freezeLimiter, &flag, 1);
      // global parameters for monitoring the simulation
      CFuint iter = Framework::SubSystemStatusStack::getActive()->getNbIter();
      CFreal res  = Framework::SubSystemStatusStack::getActive()->getResidual();    
      CudaEnv::copyHost2Dev(&dco->currIter, &iter, 1);
      CudaEnv::copyHost2Dev(&dco->currRes, &res, 1);
    }   
    
    /// copy the local configuration options to the device
    void copyConfigOptions(DeviceConfigOptions<NOTYPE>* dco) 
    {
      // consider to copy to constant memory
      dco->limitIter = _limitIter;
      dco->limitRes = _limitRes;
      dco->gradientFactor = _gradientFactor;
      dco->freezeLimiter = static_cast<CFuint>(_freezeLimiter);
      // global parameters for monitoring the simulation
      dco->currIter = Framework::SubSystemStatusStack::getActive()->getNbIter();
      dco->currRes = Framework::SubSystemStatusStack::getActive()->getResidual();
    }   
  #endif
  
    /**
     * Constructor
     */
    LeastSquareP1PolyRec3D(const std::string& name);
  
    /**
     * Default destructor
     */
    virtual ~LeastSquareP1PolyRec3D();
  
    /**
     * Configure the object
     */
    virtual void configure ( Config::ConfigArgs& args );
  
    /**
     * Returns the DataSocket's that this numerical strategy needs as sinks
     * @return a vector of SafePtr with the DataSockets
     */
    virtual std::vector<Common::SafePtr<Framework::BaseDataSocketSink> >
    needsSockets();
  
    /**
     * Compute the gradients
     */
    virtual void computeGradients();
  
    /**
     * Set up the private data
     */
    virtual void setup();
  
    /**
     * Update the weights when nodes are moving
     */
    virtual void updateWeights();
  
  protected:
  
    /**
     * Extrapolate the solution in the face quadrature points
     */
    virtual void extrapolateImpl(Framework::GeometricEntity* const face);
  
    /**
     * Extrapolate the solution in the face quadrature points
     */
    virtual void extrapolateImpl(Framework::GeometricEntity* const face,
                                 CFuint iVar, CFuint leftOrRight);
  
  protected:
  
    /// socket for stencil
    Framework::DataSocketSink<std::vector<Framework::State*> > socket_stencil;
  
    /// socket for weights
    Framework::DataSocketSink<CFreal> socket_weights;
  
    /// socket for uX values
    Framework::DataSocketSink<CFreal> socket_uX;
  
    /// socket for uY values
    Framework::DataSocketSink<CFreal> socket_uY;
  
    /// socket for uZ values
    Framework::DataSocketSink<CFreal> socket_uZ;
  
    RealVector  _l11;
  
    RealVector  _l12;
  
    RealVector  _l13;
  
    RealVector  _l22;
  
    RealVector  _l23;
  
    RealVector  _l33;
  
    RealVector  _lf1;
  
    RealVector  _lf2;
  
    RealVector  _lf3;
  
  }; // end of class LeastSquareP1PolyRec3D
  
  //////////////////////////////////////////////////////////////////////////////
  
  #ifdef CF_HAVE_CUDA
  
  template <typename PHYS>
  void LeastSquareP1PolyRec3D::DeviceFunc<PHYS>::computeGradients
  (const CFreal* state,  const CFreal* node, KernelData<CFreal> *const kd, CellData::Itr* cell)
  {
    CFreal l11 = 0.0;
    CFreal l12 = 0.0;
    CFreal l22 = 0.0;
    CFreal l13 = 0.0;
    CFreal l23 = 0.0;
    CFreal l33 = 0.0;
    CFreal lf1[PHYS::NBEQS];
    CFreal lf2[PHYS::NBEQS];
    CFreal lf3[PHYS::NBEQS];
    for (CFuint i = 0; i< PHYS::NBEQS; ++i) {
      lf1[i] = 0.;
      lf2[i] = 0.;
      lf3[i] = 0.;
    }  
    
    const CFreal* states = &kd->states[0];
    const CFreal* ghostStates = &kd->ghostStates[0];
    const CFreal* centerNodes = &kd->centerNodes[0];
    const CFreal* ghostNodes = &kd->ghostNodes[0];
    
    // loop over the neighbor cells belonging to the chosen stencil
    const CFuint stencilSize = cell->getStencilSize();
    for(CFuint in = 0; in < stencilSize; ++in) {
      const CFuint cellID = cell->getNeighborID(in);
      const CFint stype = cell->getNeighborType(in);
      // here a check on stype != 0 should not be needed since all
      // partition faces are neglected when stencil is constructed
      const CFreal* nodeLast = (stype > 0) ? &centerNodes[cellID*PHYS::DIM] : &ghostNodes[cellID*PHYS::DIM];  
      const CFreal deltaR = getDistance<PHYS::DIM, CFreal, CFreal>(node, nodeLast);
      const CFreal weight = 1.0/deltaR;
      const CFreal dx = weight*(nodeLast[XX] - node[XX]);
      const CFreal dy = weight*(nodeLast[YY] - node[YY]);
      const CFreal dz = weight*(nodeLast[ZZ] - node[ZZ]);
      l11 += dx*dx;
      l12 += dx*dy;
      l22 += dy*dy;
      l13 += dx*dz;
      l23 += dy*dz;
      l33 += dz*dz;
      
      const CFreal* neighState = (stype > 0) ? &states[cellID*PHYS::NBEQS] : &ghostStates[cellID*PHYS::NBEQS];  
      for (CFuint iVar = 0; iVar < PHYS::NBEQS; ++iVar) {
        const CFreal du = weight*(neighState[iVar] - state[iVar]);
        lf1[iVar] += dx*du;
        lf2[iVar] += dy*du;
        lf3[iVar] += dz*du;
      }
    }
    
    const CFreal det = l11*l22*l33 - l11*l23*l23 - l12*l12*l33
      + l12*l13*l23 + l13*l12*l23 - l13*l13*l22;
    
    const CFreal linv11 =   l22*l33 - l23*l23;
    const CFreal linv22 =   l11*l33 - l13*l13;
    const CFreal linv33 =   l11*l22 - l12*l12;
    const CFreal linv12 = -(l12*l33 - l13*l23);
    const CFreal linv13 =   l12*l23 - l13*l22;
    const CFreal linv23 = -(l11*l23 - l13*l12);
  
    CFreal* uX = &kd->uX[0];
    CFreal* uY = &kd->uY[0];
    CFreal* uZ = &kd->uZ[0];
  
    // A cure to the singularites in calculating the determinant
    const CFuint starts = cell->getCellID()*PHYS::NBEQS;
    if (abs(det) > 1e-16) { // maybe 1e-12 would be more conservative...
      const CFreal invDet = 1./det;
      for (CFuint i = 0; i < PHYS::NBEQS; ++i) {
        const CFuint gradx = starts + i;
        uX[gradx] = (linv11*lf1[i] + linv12*lf2[i] + linv13*lf3[i])*invDet;
        uY[gradx] = (linv12*lf1[i] + linv22*lf2[i] + linv23*lf3[i])*invDet;
        uZ[gradx] = (linv13*lf1[i] + linv23*lf2[i] + linv33*lf3[i])*invDet;
      }
    }
    else {
      for (CFuint i = 0; i < PHYS::NBEQS; ++i) {
        const CFuint gradx = starts + i;
        uX[gradx] = uY[gradx] = uZ[gradx] = 0.;
      }
    }
  }
  
  //////////////////////////////////////////////////////////////////////////////
  
  template <typename PHYS>
  void LeastSquareP1PolyRec3D::DeviceFunc<PHYS>::extrapolateOnFace
  (FluxData<PHYS>* currFd, const CFreal* coord, const CFreal* uX, 
   const CFreal* uY, const CFreal* uZ, const CFreal* limiter)
  {
    // gradient is not computed for ghost states (the corresponding inner gradient is taken instead)
    const CFuint stateIDL = currFd->getStateID(LEFT);
    const CFuint stateIDR = (!currFd->isBFace()) ? currFd->getStateID(RIGHT) : stateIDL;
    const CFreal* stateR = currFd->getState(RIGHT);
    const CFreal* stateL = currFd->getState(LEFT);
    CFreal* rstateR = currFd->getRstate(RIGHT);
    CFreal* rstateL = currFd->getRstate(LEFT);
    CFreal* nodeR = currFd->getNode(RIGHT);
    CFreal* nodeL = currFd->getNode(LEFT);
    
    // L/R cell coordinates are overwritten with the quadrature point 
    CFreal* rnodeR = currFd->getRnode(RIGHT);
    CFreal* rnodeL = currFd->getRnode(LEFT);
    for (CFuint i = 0; i < PHYS::DIM; ++i) {
      rnodeR[i] = rnodeL[i] = coord[i];
    }
    
    // AL: the following is more consistent in first order
    // for (CFuint i = 0; i < PHYS::DIM; ++i) {
    //   rnodeR[i] = (m_dco->gradientFactor > 0.) ? coord[i] : nodeR[i];
    //   rnodeL[i] = (m_dco->gradientFactor > 0.) ? coord[i] : nodeL[i];
    // }
    
    const CFuint startR = stateIDR*PHYS::NBEQS;
    const CFuint startL = stateIDL*PHYS::NBEQS;
    for (CFuint iVar = 0; iVar < PHYS::NBEQS; ++iVar) {
      const CFuint sR = startR + iVar;
      const CFreal duxR = uX[sR]*(coord[XX] - nodeR[XX]);
      const CFreal duyR = uY[sR]*(coord[YY] - nodeR[YY]);
      const CFreal duzR = uZ[sR]*(coord[ZZ] - nodeR[ZZ]);
      const CFreal duTotR = duxR + duyR + duzR;
      rstateR[iVar] = stateR[iVar];
      rstateR[iVar] += m_dco->gradientFactor*limiter[sR]*duTotR;
      const CFuint sL = startL + iVar;
      const CFreal duxL = uX[sL]*(coord[XX] - nodeL[XX]);
      const CFreal duyL = uY[sL]*(coord[YY] - nodeL[YY]);
      const CFreal duzL = uZ[sL]*(coord[ZZ] - nodeL[ZZ]);
      const CFreal duTotL = duxL + duyL + duzL;
      rstateL[iVar] = stateL[iVar];
      rstateL[iVar] += m_dco->gradientFactor*limiter[sL]*duTotL;
    }
  }
  
  //////////////////////////////////////////////////////////////////////////////
  
  template <typename PHYS>
  inline void LeastSquareP1PolyRec3D::DeviceFunc<PHYS>::extrapolateOnFace
  (const CFuint iVar, FluxData<PHYS>* currFd, const CFreal* coord, const CFreal* uX, 
   const CFreal* uY, const CFreal* uZ, const CFreal* limiter)
  {
    // gradient is not computed for ghost states (the corresponding inner gradient is taken instead)
    const CFuint sL   = currFd->getStateID(LEFT)*PHYS::NBEQS + iVar;
    const CFreal* nodeL = currFd->getNode(LEFT);
    const CFreal duxL = uX[sL]*(coord[XX] - nodeL[XX]);
    const CFreal duyL = uY[sL]*(coord[YY] - nodeL[YY]);
    const CFreal duzL = uZ[sL]*(coord[ZZ] - nodeL[ZZ]);
    const CFreal duTotL = duxL + duyL + duzL;
    currFd->getRstate(LEFT)[iVar] = currFd->getState(LEFT)[iVar];
    currFd->getRstate(LEFT)[iVar] += m_dco->gradientFactor*limiter[sL]*duTotL; 
  }
  
  #endif
        
  //////////////////////////////////////////////////////////////////////////////
        
      } // namespace FiniteVolume
  
    } // namespace Numerics
  
  } // namespace COOLFluiD
  
  //////////////////////////////////////////////////////////////////////////////
  
  #endif // COOLFluiD_Numerics_FiniteVolume_LeastSquareP1PolyRec3D_hh
  