CutCellMomentsImplem.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

#ifndef _CUTCELLMOMENTSIMPLEM_H_
#define _CUTCELLMOMENTSIMPLEM_H_


#if defined(CH_Darwin) && defined(__GNUC__) && ( __GNUC__ == 3 )
// deal with the broken isnan()/isinf() in GCC on MacOS
#include <unistd.h>
#define _GLIBCPP_USE_C99 1
#endif

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>

#include "MayDay.H"

#include "NoRefinement.H"
#include "LSProblem.H"

#include "NamespaceHeader.H"

// Null constructor
template <int dim> inline CutCellMoments<dim>::CutCellMoments()
{
}

inline CutCellMoments<1>::CutCellMoments()
{
}

// Copy constructor
template <int dim> inline CutCellMoments<dim>::CutCellMoments(const CutCellMoments<dim>& a_cutCellMoments)
  :m_moments(a_cutCellMoments.m_moments),
   m_EBmoments(a_cutCellMoments.m_EBmoments),
   m_bdCutCellMoments(a_cutCellMoments.m_bdCutCellMoments),
   m_iFData(a_cutCellMoments.m_iFData),
   m_bdCCOn(a_cutCellMoments.m_bdCCOn),
   m_residual(a_cutCellMoments.m_residual)
{
}

// Copy constructor
inline CutCellMoments<1>::CutCellMoments(const CutCellMoments<1>& a_thisCutCellMoments)

  :m_moments(a_thisCutCellMoments.m_moments),
   m_iFData(a_thisCutCellMoments.m_iFData)
{
}

// Use this for initializing
template <int dim> inline CutCellMoments<dim>::CutCellMoments(const IFData<dim>& a_info)
  :m_iFData(a_info)
{
  CH_TIME("CutCellMoments::constructor");

  m_bdCCOn = false;

  for (int hilo = 0; hilo < 2; ++hilo)
  {
    for (int idir = 0; idir < dim; ++idir)
    {
      Iv2 bdId;
      bdId[BDID_DIR] = idir;
      bdId[BDID_HILO] = hilo;

      IFData<dim-1> reducedInfo(a_info,idir,hilo);
      CutCellMoments<dim-1>bdCutCellMoments(reducedInfo);

      m_bdCutCellMoments[bdId] = bdCutCellMoments;

      // Notice whether at least one lower dimensional cutCell is on the
      // interface
      if (reducedInfo.m_allVerticesOn)
      {
        m_bdCCOn = true;
      }
    }
  }
}

// This constructor is used in the recursion
inline CutCellMoments<1>::CutCellMoments(const IFData<1> & a_info)
  :m_iFData(a_info)
{
}

// Destructor
template <int dim> inline CutCellMoments<dim>::~CutCellMoments()
{
}

inline CutCellMoments<1>::~CutCellMoments()
{
}

template <int dim> inline const CutCellMoments<dim - 1> CutCellMoments<dim>::getBdCutCellMoments(const Iv2 & a_bdId) const
{
  typename BdCutCellMoments::const_iterator it = m_bdCutCellMoments.find(a_bdId);

  if (it == m_bdCutCellMoments.end())
  {
    MayDay::Abort("Can't find this bdId in m_bdCutCellMoments");
  }

  return it->second;
}

// Solve
template <int dim> inline void CutCellMoments<dim>::computeMoments(const int           & a_order,
                                                                   const int           & a_degreeP,
                                                                   const bool          & a_useConstraints,
                                                                   RefinementCriterion & a_refinementCriterion)
{
  CH_TIME("CutCellMoments::ComputeMoments");

   Vector<Real> RNorm(3);
   for (int i = 0; i < 3; i++)
   {
     RNorm[i] = LARGEREALVAL;
   }

   for (int i = 0; i < a_degreeP + 1; i++)
   {
     m_residual.push_back(RNorm);
   }

   if (m_iFData.m_allVerticesOut)
   {
     for (int iDegree = a_degreeP; iDegree >= 0; --iDegree)
     {
       LSProblem<dim> lsp(iDegree,a_useConstraints);

       // Fill moments
       const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
       for (typename PthMomentLoc::const_iterator it=monDegreeP.begin(); it != monDegreeP.end(); ++it)
       {
         m_EBmoments[it->first] = 0.0;
       }

       const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
       for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
       {
         m_moments[it->first] = 0.0;
       }
     }
   }
   else if (m_iFData.m_allVerticesIn && !m_bdCCOn)
   {
     for (int iDegree = a_degreeP; iDegree >= 0; --iDegree)
     {
       LSProblem<dim> lsp(iDegree,a_useConstraints);

       // Fill moments of degree P and P-1
       const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
       for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
       {
         m_EBmoments[it->first] = 0.0;
       }

       const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
       for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
       {
         m_moments[it->first] = fullCellQuadrature(it->first);
       }
     }
   }
   else
   {
     for (typename BdCutCellMoments::iterator it = m_bdCutCellMoments.begin(); it != m_bdCutCellMoments.end(); ++it)
     {
       NoRefinement noRefinement;
       it->second.computeMoments(a_order,a_degreeP+1,a_useConstraints,noRefinement);
     }

     for (int iDegree = a_degreeP; iDegree >= 0; --iDegree)
     {
       // Make a LSProb
       LSProblem<dim> lsp(a_order,iDegree,a_useConstraints,m_iFData.m_normal);

       Vector<Real> rhs = computeRhs(lsp,a_order);

       if (a_useConstraints)
       {
         lsp.computeBounds(m_iFData.m_dx,*this);
       }

       // Solve the problem and return residual
       lsp.invertNormalEq(rhs,m_residual[iDegree]);

       // Fill moments
       const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
       for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
       {
         m_EBmoments[it->first] = lsp.getUnknown(it->second);
       }

       const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
       for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
       {
         m_moments[it->first] = lsp.getUnknown(it->second
                              + lsp.getNumberDegP());
       }
     }
   }

   Vector<int>  refineInDir(dim);
   Vector<Real> vectDx     (dim);

   for (int idir = 0; idir < dim; idir++)
   {
     vectDx[idir] = m_iFData.m_dx[idir];
   }

   if (a_refinementCriterion.doRefine(refineInDir,dim,vectDx,m_residual))
   {
     // Refine creates a vector of 2^dim cutcellmoments corresponding to
     // the refined cells inside the current cell
     Vector<CutCellMoments<dim> > refinedCCMoms = refine(a_order,a_degreeP,a_useConstraints,a_refinementCriterion,refineInDir);

     // addMomentMaps add up the refined moments in the current
     // CutCellMoments maps
     addMomentMaps(refinedCCMoms,a_degreeP,a_useConstraints);

     // computeResiduals(a_order,a_degreeP,a_useConstraints);
     computeResiduals(refinedCCMoms,a_degreeP);
   }
}

// Integrate along line segments aligned in a coordinalte direction
inline void CutCellMoments<1>::computeMoments(const int           & a_order,
                                              const int           & a_degree,
                                              const bool          & a_useConstraints,
                                              RefinementCriterion & a_refinementCriterion)
{
  int lo = 0;
  int hi = 1;
  int loSign = m_iFData.m_cornerSigns[lo];
  int hiSign = m_iFData.m_cornerSigns[hi];

  // If entire edge out of the fluid, then moments = 0.0
  if (loSign <= ON && hiSign<= ON)
  {
    for (int iDegree = 0; iDegree <= a_degree; ++iDegree)
    {
      // Definition of m_moments typedef requires that we define degree as
      // a oneTuple:
      IndexTM<int,1> degree;
      degree[0] = iDegree;
      m_moments[degree] = 0.0;
    }
  }
  else
  {
    // Assign loPt and hiPt with respect to coordinates centered at cell-
    // center
    Real loPt = LARGEREALVAL;
    Real hiPt = LARGEREALVAL;

    if (loSign >= ON)
    {
      loPt = -0.5 * m_iFData.m_dx;
    }
    else
    {
      loPt =  0.5 * m_iFData.m_intersection * m_iFData.m_dx;
    }

    if (hiSign >= ON)
    {
      hiPt =  0.5 * m_iFData.m_dx;
    }
    else
    {
      hiPt =  0.5 * m_iFData.m_intersection * m_iFData.m_dx;
    }

    CH_assert(-0.5*m_iFData.m_dx <= loPt && loPt <= 0.5*m_iFData.m_dx);
    CH_assert(-0.5*m_iFData.m_dx <= hiPt && hiPt <= 0.5*m_iFData.m_dx);

    // Integrate x^degree over the line segment[loPt,hiPt]
    computeMomentsUsingBinomial(loPt,hiPt,loSign,hiSign,a_degree);
  }
}

inline void CutCellMoments<1>::simpleComputeMoments(const Real & a_loPt,
                                                    const Real & a_hiPt,
                                                    const int  & a_degree)
{
  for (int iDegree = 0; iDegree <= a_degree; ++iDegree)
  {
    //definition of m_moments typedef requires that we define degree thus:
    IndexTM<int,1>degree;
    degree[0] = iDegree;
    m_moments[degree] = pow(a_hiPt,iDegree + 1) - pow(a_loPt,iDegree +1);
    m_moments[degree] /= (iDegree + 1);
  }
}

inline void CutCellMoments<1>::computeMomentsUsingBinomial(const Real & a_loPt,
                                                           const Real & a_hiPt,
                                                           const int  & a_loSign,
                                                           const int  & a_hiSign,
                                                           const int  & a_degree)
{
  for (int iDegree = 0; iDegree <= a_degree; ++iDegree)
  {
    Real epsilon = a_hiPt - a_loPt;

    // Definition of m_moments typedef requires that we define degree thus:
    IndexTM<int,1> degree;
    degree[0] = iDegree;
    m_moments[degree] = 0.0;

    // New method without substracting higher order terms
    for (int j = 1; j <= iDegree + 1; j++)
    {
      int bigger = j;
      int smaller = j;
      if (iDegree + 1 - j > j)
      {
        bigger = iDegree + 1 - j;
      }
      else
      {
        smaller = iDegree + 1 - j;
      }

      int numerator = 1;
      for (int i = bigger + 1; i <= iDegree + 1; ++i)
      {
        numerator *= i;
      }

      int denominator = 1;
      for (int i = 1; i <= smaller; ++i)
      {
        denominator *= i;
      }

      Real factor = numerator / denominator;
      if (a_loSign >= ON)
      {
        m_moments[degree] += factor * pow(a_loPt,iDegree + 1 - j) * pow(epsilon,j);
      }
      else if (a_hiSign >= ON)
      {
        m_moments[degree] -= factor * pow(a_hiPt,iDegree + 1 - j) * pow(epsilon,j) * pow(-1.0,j);
      }
    }

    m_moments[degree] /= iDegree + 1;
  }
}

template <int dim> inline Real CutCellMoments<dim>::fullCellQuadrature(const IvDim & a_mono)
{
  Real moment = 1.0;
  RvDim vectDx = m_iFData.m_dx;

  for (int idir = 0; idir < dim; ++idir)
  {
    Real idirInt = pow( 0.5*vectDx[idir],a_mono[idir] + 1)
                 - pow(-0.5*vectDx[idir],a_mono[idir] + 1);

    idirInt /= (a_mono[idir] + 1);
    moment *= idirInt;
  }

  return moment;
}

template <int dim> inline Vector<Real> CutCellMoments<dim>::computeRhs(LSProblem<dim> & a_lsp,
                                                                       const int      & a_order)
{
  // Resize rhs
  int numEq = dim*a_lsp.getNumberDegP();
  Vector<Real> rhs(numEq);

  Vector< IndexTM<Real,dim> > gradNorm = m_iFData.m_gradNormal;

  // For each moment iterate thru bd CutCellMoments incrementing same comp
  // of rhs
  const LocPthMoment& locMap = a_lsp.getLocMonomialMapDegreeP();
  for (typename LocPthMoment::const_iterator it = locMap.begin(); it != locMap.end(); ++it)
  {
    int jth = it->first;
    IvDim mono = it->second;

    int hiSide = 1;
    int loSide = 0;
    Iv2 bdId;

    for (int idir = 0; idir < dim; ++idir)
    {
      // Which lower dimensional monomial corresponds (mono,j)
      IndexTM<int,dim-1> mono1Less;
      for (int jdir = 0; jdir < dim; ++jdir)
      {
        if (jdir < idir)
        {
          mono1Less[jdir] = mono[jdir];
        }
        else if (jdir > idir)
        {
          mono1Less[jdir-1] = mono[jdir];
        }
      }

      bdId[0] = idir;
      bdId[1] = hiSide;

      Real hiMom = m_bdCutCellMoments[bdId].m_moments[mono1Less];

      bdId[1] = loSide;

      Real loMom = m_bdCutCellMoments[bdId].m_moments[mono1Less];
      int exponent = it->second[idir];

      Real hiFactor = pow(0.5*m_iFData.m_dx[idir],exponent);
      Real loFactor = pow(-0.5*m_iFData.m_dx[idir],exponent);

      rhs[(dim*jth) + idir] = hiMom*hiFactor - loMom*loFactor;

      // Add the gradNormal term
      if (a_order < 0 || a_order > 1)
      {
        MayDay::Abort("CutCellMoments::computeMoments not set up for orders not 0 or 1");
      }

      if (a_order == 1)
      {
        // Add the gradNormal term
        for (int jdir = 0; jdir < dim; jdir++)
        {
          IvDim monoPlus1 = mono;
          monoPlus1[jdir] += 1;

          if (m_EBmoments.find(monoPlus1) != m_EBmoments.end())
          {
            Real moment = m_EBmoments[monoPlus1];
            rhs[(dim*jth) + idir] += gradNorm[idir][jdir] * moment;
          }
        }
      }
    }
  }

  return rhs;
}

template <int dim> inline Vector< CutCellMoments<dim> > CutCellMoments<dim>::refine(const int           & a_order,
                                                                                    const int           & a_degreeP,
                                                                                    const bool          & a_useConstraints,
                                                                                    RefinementCriterion & a_refinementCriterion,
                                                                                    const Vector<int>   & a_refineInDir)
{
  CH_assert(m_iFData.m_changingDir.size() == dim);

  Vector<CutCellMoments<dim> > refinedCutCellMoments;

  // This function refines the current cell by a factor of two in all
  // directions where a_refineInDir is non-zero
  IndexTM<Real,dim> refinedDx = m_iFData.m_dx;

  // Current cell iterator and maximum iterator
  IndexTM<int,dim> curIter;
  IndexTM<int,dim> maxIter;

  for (int idir = 0; idir < dim; idir++)
  {
    if (a_refineInDir[idir] != 0)
    {
      // This direction will be divided in two
      refinedDx[idir] /= 2.0;
      maxIter[idir]    = 1;
    }
    else
    {
      // This direction is undivided
      maxIter[idir]    = 0;
    }

    // Initial the start to zero
    curIter[idir] = 0;
  }

  // Iterate through all the subcells
  while (1)
  {
    // Computes the cell center for the iCell refined cell
    IndexTM<Real,GLOBALDIM> refinedCellCenter = m_iFData.m_globalCellCenter;

    for (int idir = 0; idir < dim; idir++)
    {
      int globalDir = m_iFData.m_changingDir[idir];

      if (a_refineInDir[idir] != 0)
      {
        int sign = 2*curIter[idir] - 1;
        refinedCellCenter[globalDir] += sign * 0.5*refinedDx[idir];
      }
    }

    // Creates the IF Data for the refined cell, it uses the global
    // constructor of IFData and recomputes all the intersection points
    // and normal/grad normal from the implicit function
    IFData<dim> refinedIFData(*m_iFData.m_function,refinedDx,refinedCellCenter,m_iFData.m_normalDir);

    // Creates a cutcellmoments for the refined cell
    CutCellMoments<dim> refinedCCMom(refinedIFData);

    // Compute the moments for the refined cell
    refinedCCMom.computeMoments(a_order,a_degreeP,a_useConstraints,a_refinementCriterion);
    refinedCutCellMoments.push_back(refinedCCMom);

    // Move to the next subcell (if there is one)
    int idir;
    for (idir = 0; idir < dim; idir++)
    {
      if (curIter[idir] < maxIter[idir])
      {
        // If the current index is in range, increment it and continue
        curIter[idir]++;
        break;
      }
      else
      {
        // If the current index at the maximum, reset it and move to the next
        // index
        curIter[idir] = 0;
      }
    }

    // All the subcells have been done
    if (idir == dim)
    {
      break;
    }
  }

  return refinedCutCellMoments;
}

template <int dim> inline void CutCellMoments<dim>::addMomentMaps(const Vector<CutCellMoments<dim> > & a_refinedCutCellVector,
                                                                  const int                          & a_degreeP,
                                                                  const bool                         & a_useConstraints)
{
  int numberOfRefinedCC = a_refinedCutCellVector.size();
  CH_assert(numberOfRefinedCC>0);

  IndexTM<int,2> bdId;

  // Pick one refined CutCellMoments and use its maps to initialize the maps
  // of the current CutCellMoment to zero the chosen cutcellmoment need to
  // have initialized boundary moments, i.e., be irregular
  CutCellMoments<dim> cutCell;
  int k = 0;

  while ((a_refinedCutCellVector[k].m_iFData.m_allVerticesOut || a_refinedCutCellVector[k].m_iFData.m_allVerticesIn)
      && k < numberOfRefinedCC)
  {
    k++;
  }

  if (k >= numberOfRefinedCC)
  {
    MayDay::Abort("Exceeded the number of refined CutCellMoments in vector");
  }
  else
  {
    cutCell = a_refinedCutCellVector[k];
  }

  initialize(cutCell);

  // Add up the moments on the refined cells to obtain the moment in the
  // coarse cell moments need to be add up in the same coordinate system
  // which is the one where the origin is the cell center of the current cell
  for (int i = 0; i < numberOfRefinedCC; i++)
  {
    // The difference between the refined cell center and the current cell
    // center is needed to change the coordinate system of the refined moments
    // to the coordinate system where they are needed.
    CutCellMoments<dim> refinedCutCell = a_refinedCutCellVector[i];
    IndexTM<int,dim> hilo = LARGEINTVAL*IndexTM<int,dim>::Unit;
    IndexTM<Real,dim> refinedCenterDelta = IndexTM<Real,dim>::Zero;

    int ii = i;
    for (int j = 0; j < dim; ++j)
    {
      hilo[j] = (ii & 1);
      ii = ii >> 1;
    }

    IndexTM<int,dim> sign = 2*hilo-1;

    for (int idir = 0; idir < dim; idir++)
    {
      refinedCenterDelta[idir] = 0.25*sign[idir]*m_iFData.m_dx[idir];
    }

    // Add up the moments in the volume
    addMoments(m_moments,refinedCutCell.m_moments,refinedCenterDelta);
    addMoments(m_EBmoments,refinedCutCell.m_EBmoments,refinedCenterDelta);

    // Add up the moments on the boundary
    for (int idir = 0; idir < dim; idir++)
    {
      bdId[0] = idir;

      // The difference in cell centers in dim-1 is needed to change the
      // coordinate system of the refined moments to the coordinate
      // system where they are needed.
      IndexTM<Real,dim-1> localRefinedCenterDelta;
      IndexTM<int,dim-1> localHilo;

      for (int jdir = 0; jdir < dim; jdir++)
      {
        if (jdir < idir)
        {
          localRefinedCenterDelta[jdir] = refinedCenterDelta[jdir];
          localHilo[jdir] = hilo[jdir];
        }
        else if (jdir > idir)
        {
          localRefinedCenterDelta[jdir-1] = refinedCenterDelta[jdir];
          localHilo[jdir-1] = hilo[jdir];
        }
      }

      if (hilo[idir] != LARGEINTVAL)
      {
        bdId[1] = hilo[idir];

        // Drop one dimension to add up the boundary moments dropping
        // one dimension allows an easier implementation as 1D is a
        // particular case
        (refinedCutCell.m_bdCutCellMoments[bdId]).addBdMoments(
            m_bdCutCellMoments[bdId],
            refinedCutCell.m_iFData,
            a_degreeP,
            a_useConstraints,
            localRefinedCenterDelta,
            localHilo);
      }
      else
      {
        for (int side = 0; side < 2; side++)
        {
          bdId[1] = side;
          (refinedCutCell.m_bdCutCellMoments[bdId]).addBdMoments(
              m_bdCutCellMoments[bdId],
              refinedCutCell.m_iFData,
              a_degreeP+1,
              a_useConstraints,
              localRefinedCenterDelta,
              localHilo);
        }
      }
    }
  }
}

template <int dim> inline void CutCellMoments<dim>::addMoments(PthMoment               & a_momentMap,
                                                               PthMoment               & a_refinedMomentMap,
                                                               const IndexTM<Real,dim> & a_refinedCenterDelta)
{
  // Iterate through the REFINEDmap and add the refined moment in the right
  // coordinate system to the coarse map
  for (typename PthMoment::const_iterator it = a_refinedMomentMap.begin(); it != a_refinedMomentMap.end(); ++it)
  {
    IndexTM<int,dim> monomial = it->first;
    Real moment = changeMomentCoordinates(a_refinedMomentMap,monomial,a_refinedCenterDelta);

    a_momentMap[it->first] += moment;
  }
}

// This function is called to fill the boundary moments from the higher
// dimension: it fills the moments in a_coarseCutCell using either the
// moments in the current cutcell map or the full cell/covered cell answers
// depending on the value of a_iFData.m_allVerticesIn/Out
template <int dim> inline void CutCellMoments<dim>::addBdMoments(CutCellMoments<dim>     & a_coarseBdCutCell,
                                                                 const IFData<dim+1>     & a_iFData,
                                                                 const int               & a_degreeP,
                                                                 const bool              & a_useConstraints,
                                                                 const IndexTM<Real,dim> & a_refinedCenterDelta,
                                                                 const IndexTM<int,dim>  & a_hilo)
{
  typedef map<IvDim,int,LexLT <IvDim> > PthMomentLoc;

  if (a_iFData.m_allVerticesIn)
  {
    // The boundary map doesn't exist, we need to compute the full cell
    // quadrature for all monomial in order to use the changecoordinates
    // function.  Generate a full cell quadrature map.
    PthMoment fullCellMap;

    for (int iDegree = a_degreeP; iDegree >= 0; iDegree--)
    {
      // Creating a lsproblem generates the monomial map that is
      // needed to fill the moment maps
      LSProblem<dim> lsp(iDegree,a_useConstraints);
      const PthMomentLoc& momMap = lsp.getMonomialLocMapDegreePLess1();

      for (typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
      {
        fullCellMap[it->first] = fullCellQuadrature(it->first);
      }
    }

    for (int iDegree = a_degreeP; iDegree >= 0; iDegree--)
    {
      LSProblem<dim> lsp(iDegree,a_useConstraints);

      // Change the coordinates of the refined moments using the full
      // cell quadrature map.  Add up the obtained moment to the
      // coarse map
      const PthMomentLoc& momMap = lsp.getMonomialLocMapDegreePLess1();

      for (typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
      {
        Real bdMoment = getBdMoment(it->first,a_iFData,a_refinedCenterDelta,fullCellMap);
        a_coarseBdCutCell.m_moments[it->first] += bdMoment;
      }
    }
  }
  else
  {
    for (int iDegree = a_degreeP; iDegree >= 0; iDegree--)
    {
      LSProblem<dim> lsp(iDegree,a_useConstraints);

      // Get the moments from the refined cell map + change coordinates
      // function if the refined cell is irregular or they are 0 as all
      // vertices are out
      const PthMomentLoc& momMap = lsp.getMonomialLocMapDegreePLess1();

      for (typename PthMomentLoc::const_iterator it = momMap.begin(); it != momMap.end(); ++it)
      {
        Real bdMoment = getBdMoment(it->first,a_iFData,a_refinedCenterDelta);
        a_coarseBdCutCell.m_moments[it->first] += bdMoment;
      }
    }
  }

  for (int iDegree = a_degreeP; iDegree >= 0; iDegree--)
  {
    LSProblem<dim> lsp(iDegree,a_useConstraints);

    // EB moments are either found in the EB moment map or equals to 0
    // (allIn or allOut)
    const PthMomentLoc& momMapP = lsp.getMonomialLocMapDegreeP();

    for (typename PthMomentLoc::const_iterator it = momMapP.begin(); it != momMapP.end(); ++it)
    {
      Real bdEBMoment = getBdEBMoment(it->first,a_iFData,a_refinedCenterDelta);
      a_coarseBdCutCell.m_EBmoments[it->first] += bdEBMoment;
    }
  }

  IndexTM<int,2> bdId;

  for (int idir = 0; idir < dim; idir++)
  {
    bdId[0] = idir;

    // The refined cell center in dim-1 is needed to change the coordinate
    // system of the refined moments to the coordinate system where they are
    // needed.
    IndexTM<Real,dim-1> localRefinedCellCenter;
    IndexTM<int,dim-1> localHilo;

    for (int jdir = 0; jdir < dim; jdir++)
    {
      if (jdir < idir)
      {
        localRefinedCellCenter[jdir] = a_refinedCenterDelta[jdir];
        localHilo[jdir] = a_hilo[jdir];
      }
      else if (jdir > idir)
      {
        localRefinedCellCenter[jdir-1] = a_refinedCenterDelta[jdir];
        localHilo[jdir-1] = a_hilo[jdir];
      }
    }

    if (a_hilo[idir] != LARGEINTVAL)
    {
      bdId[1] = a_hilo[idir];
      // Drop one dimension to add up the boundary moments. Dropping one
      // dimension allows an easier implementation as 1D is a particular
      // case
      (m_bdCutCellMoments[bdId]).addBdMoments(a_coarseBdCutCell.m_bdCutCellMoments[bdId],m_iFData,a_degreeP+1,a_useConstraints,localRefinedCellCenter,localHilo);
    }
    else
    {
      for (int side = 0; side < 2; side++)
      {
        bdId[1]=side;
        (m_bdCutCellMoments[bdId]).addBdMoments(a_coarseBdCutCell.m_bdCutCellMoments[bdId],m_iFData,a_degreeP,a_useConstraints,localRefinedCellCenter,localHilo);
      }
    }
  }
}

void CutCellMoments<1>::addBdMoments(CutCellMoments<1>     & a_coarseBdCutCell,
                                     const IFData<2>       & a_iFData,
                                     const int             & a_degreeP,
                                     const bool            & a_useConstraints,
                                     const IndexTM<Real,1> & a_refinedCenterDelta,
                                     const IndexTM<int,1>  & a_localHilo)
{
  OneDMoments fullCellMap;

  if (a_iFData.m_allVerticesIn)
  {
    IndexTM<int,1> degree;

    for (int iDegree = 0; iDegree <= a_degreeP; ++iDegree)
    {
      degree[0] = iDegree;
      Real loPt = -0.5*m_iFData.m_dx;
      Real hiPt = 0.5*m_iFData.m_dx;
      fullCellMap[degree] = pow(hiPt,degree[0] + 1) - pow(loPt,degree[0] +1);
      fullCellMap[degree] /= degree[0]+1;
    }
  }

  for (int iDegree = 0; iDegree <= a_degreeP; ++iDegree)
  {
    IndexTM<int,1>degree;
    degree[0] = iDegree;
    Real bdMoment;
    if (a_iFData.m_allVerticesIn)
    {
      bdMoment = getBdMoment(degree,a_iFData,a_refinedCenterDelta,fullCellMap);
    }
    else
    {
      bdMoment = getBdMoment(degree,a_iFData,a_refinedCenterDelta);
    }

    a_coarseBdCutCell.m_moments[degree] += bdMoment;
  }
}

template <int dim> inline Real CutCellMoments<dim>::changeMomentCoordinates(PthMoment        & a_refinedMomentMap,
                                                                            const IndexTM<int,dim>  & a_monomial,
                                                                            const IndexTM<Real,dim> & a_refinedCenterDelta)
{
  // This function outputs the moment in the coordinate system where the
  // origin is the current cell center, given the moment map of the refined
  // moments computed in the coordinate system where the origin is the
  // refined cell center it uses a Taylor development around 0 of the function
  // (x+xc)^p1*(y+yc)^p2*...
  Real moment = 0.0;

  // Degree is the degree of the monomial
  int degree = 0;

  for (int i = 0; i < dim; i++)
  {
    degree += a_monomial[i];
  }

  // Loop over the possible values of r in the development of the Taylor series
  // The Taylor series has a finite expansion as the function's derivatives
  // are zero for r>degree
  for (int r = 0; r <= degree; r++)
  {
    if (r >= 0)
    {
      // Generate all the possible monomials of degree r and add up the moment
      IndexTM<int,dim> derivative;

      for (int idir = 0; idir < dim; ++idir)
      {
        derivative[idir] = 0;
      }

      while (true)
      {
        for (int j = 1; j < dim-1; ++j)
        {
          int t = r;
          for (int k = j+1; k < dim; ++k)
          {
            t -= derivative[k];
          }

          if (derivative[j] >  t)
          {
            derivative[j+1] += 1;
            derivative[j  ]  = 0;
          }
          else
          {
            break;
          }
        }

        if (derivative[dim-1] > r)
        {
          break;
        }

        derivative[0] = r;

        for (int j = 1; j < dim; ++j)
        {
          derivative[0] -= derivative[j];
        }

        // Add up the relevant term of the refined moment map to the map
        Real coeff = 1.0;
        Real factorial = 1.0;

        for (int idir = 0; idir < dim; idir++)
        {
          for (int j = 0; j < derivative[idir]; j++)
          {
            coeff *= a_monomial[idir] - j;
            factorial *= j+1;
          }

          // This is to prevent the computation of 0*0^-1 which result
          // in a nan
          if (a_monomial[idir]-derivative[idir] < 0)
          {
            if (coeff != 0)
            {
              MayDay::Abort("Coeff should be zero when the derivative has higher coefficient than the monomial");
            }
          }
          else
          {
            coeff *= pow(a_refinedCenterDelta[idir],a_monomial[idir]-derivative[idir]);
          }
        }

        moment += coeff * a_refinedMomentMap[derivative] / factorial;

        // Increments to get to the next derivative
        derivative[1] += 1;
      }
    }
  }

  return moment;
}

inline Real CutCellMoments<1>::changeMomentCoordinates(OneDMoments           & a_refinedMomentMap,
                                                       const IndexTM<int,1>  & a_monomial,
                                                       const IndexTM<Real,1> & a_refinedCenterDelta)
{
  Real moment = 0.0;
  int degree = 0;
  degree += a_monomial[0];

  // Loop over the possible values of r in the development of the Taylor series
  for (int r = 0; r <= degree; r++)
  {
    // Generate all the possible monomials of degree r and add up the moment
    IndexTM<int,1> derivative;
    derivative[0] = r;

    // Add up the relevant term of the refined moment map to the map
    Real coeff = 1.0;
    Real factorial = 1.0;

    for (int j = 0; j < derivative[0]; j++)
    {
      coeff *= a_monomial[0] - j;
      factorial *= j+1;
    }

    coeff *= pow(a_refinedCenterDelta[0],a_monomial[0]-derivative[0]);

    moment += coeff * a_refinedMomentMap[derivative] / factorial;
  }

  return moment;
}

template <int dim> inline void CutCellMoments<dim>::initialize(CutCellMoments<dim> & a_refinedCutCell)
{
  // This initialize the maps of the current cut cell moment to zero by
  // iterating through the refined maps
  initializeMap(m_EBmoments,a_refinedCutCell.m_EBmoments);
  initializeMap(m_moments,a_refinedCutCell.m_moments);

  // Initialize maps on the boundary
  IndexTM<int,2> bdId;

  for (int idir = 0; idir < dim; idir++)
  {
    bdId[0] = idir;

    for (int hilo = 0; hilo < 2; hilo++)
    {
      bdId[1] = hilo;
      CutCellMoments<dim-1> refinedBdCutCell = a_refinedCutCell.getBdCutCellMoments(bdId);

      m_bdCutCellMoments[bdId].initialize(refinedBdCutCell);
    }
  }
}

inline void CutCellMoments<1>::initialize(CutCellMoments<1> & a_refinedCutCell)
{
  initializeMap(m_EBmoments,a_refinedCutCell.m_EBmoments);
  initializeMap(m_moments,a_refinedCutCell.m_moments);
}

template <int dim> inline void CutCellMoments<dim>::initializeMap(PthMomentLesserDimension & a_map1,
                                                                  PthMomentLesserDimension & a_map2)
{
  for (typename PthMomentLesserDimension::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
  {
    a_map1[it->first] = 0.0;
  }
}

inline void CutCellMoments<1>::initializeMap(OneDMoments & a_map1,
                                             OneDMoments & a_map2)
{
  for (OneDMoments::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
  {
    a_map1[it->first] = 0.0;
  }
}

template <int dim> inline void CutCellMoments<dim>::initializeMap(PthMoment & a_map1,
                                                                  PthMoment & a_map2)
{
  for (typename PthMoment::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
  {
    a_map1[it->first] = 0.0;
  }
}

template <int dim> inline void CutCellMoments<dim>::initializeMap(OneDMoments & a_map1,
                                                                  OneDMoments & a_map2)
{
  for (typename OneDMoments::const_iterator it = a_map2.begin(); it != a_map2.end(); ++it)
  {
    a_map1[it->first] = 0.0;
  }
}

template <int dim> inline Real CutCellMoments<dim>::getBdMoment(const IvDim             & a_mono,
                                                                const IFData<dim+1>     & a_iFData,
                                                                const IndexTM<Real,dim> & a_refinedCenterDelta,
                                                                PthMoment                 a_fullCellMap)
{
  // This computes the moments on the fly in the case the refined cell (in
  // dim+1) was all in or all out (in that case the boundary moments were
  // never computed) or looks for them in the moment map
  Real moment = LARGEREALVAL;

  if (a_iFData.m_allVerticesOut)
  {
    moment = 0.0;
  }
  else if (a_iFData.m_allVerticesIn)
  {
    moment = changeMomentCoordinates(a_fullCellMap,a_mono,a_refinedCenterDelta);
  }
  else
  {
    moment = changeMomentCoordinates(m_moments,a_mono,a_refinedCenterDelta);
  }

  return moment;
}

inline Real CutCellMoments<1>::getBdMoment(const IndexTM<int,1>  & a_mono,
                                           const IFData<2>       & a_iFData,
                                           const IndexTM<Real,1> & a_refinedCenterDelta,
                                           OneDMoments             a_fullCellMap)
{
  Real moment = LARGEREALVAL;

  if (a_iFData.m_allVerticesOut)
  {
    moment = 0.0;
  }
  else if (a_iFData.m_allVerticesIn)
  {
    Real loPt = -0.5*m_iFData.m_dx;
    Real hiPt =  0.5*m_iFData.m_dx;

    moment = pow(hiPt,a_mono[0] + 1) - pow(loPt,a_mono[0] +1);
    moment = changeMomentCoordinates(a_fullCellMap,a_mono,a_refinedCenterDelta);
  }
  else
  {
    moment = changeMomentCoordinates(m_moments,a_mono,a_refinedCenterDelta);
  }

  return moment;
}

template <int dim> inline Real CutCellMoments<dim>::getBdEBMoment(const IvDim             & a_mono,
                                                                  const IFData<dim+1>     & a_iFData,
                                                                  const IndexTM<Real,dim> & a_refinedCenterDelta)
{
  // EB moments are 0 if the refined cell (in dim+1) is all in or all out,
  // otherwise they are in the current cell EB moment maps
  Real EBmoment = LARGEREALVAL;

  if (a_iFData.m_allVerticesOut || a_iFData.m_allVerticesIn)
  {
    EBmoment = 0.0;
  }
  else
  {
    EBmoment = changeMomentCoordinates(m_EBmoments,a_mono,a_refinedCenterDelta);
  }

  return EBmoment;
}

inline Real CutCellMoments<1>::getBdEBMoment(const IndexTM<int,1>  & a_mono,
                                             const IFData<2>       & a_iFData,
                                             const IndexTM<Real,1> & a_refinedCenterDelta)
{
  Real EBmoment = 0.0;

  return EBmoment;
}

template <int dim> inline void CutCellMoments<dim>::computeResiduals(const int  & a_order,
                                                                     const int  & a_degreeP,
                                                                     const bool & a_useConstraints)
{
  for (int iDegree = a_degreeP; iDegree >= 0; iDegree--)
  {
    // Initialize residuals to 0
    int nNorm = 3;
    m_residual[iDegree].resize(nNorm);

    for (int i = 0; i < nNorm; i++)
    {
      m_residual[iDegree][i] = 0.0;
    }

    // Create a lsp problem
    LSProblem<dim> lsp(a_order,iDegree,a_useConstraints,m_iFData.m_normal);

    // Compute the right hand side
    Vector<Real> rhs = computeRhs(lsp,a_order);

    // Get the unknowns from the moment maps
    Vector<Real> unknowns(lsp.m_numP+lsp.m_numPLess1);

    for (typename PthMomentLoc::const_iterator it = lsp.m_monoLocP.begin(); it != lsp.m_monoLocP.end(); ++it)
    {
      unknowns[it->second] = m_EBmoments[it->first];
    }

    for (typename PthMomentLoc::const_iterator it = lsp.m_monoLocPLess1.begin(); it != lsp.m_monoLocPLess1.end(); ++it)
    {
      unknowns[it->second + lsp.m_numP] = m_moments[it->first];
    }

    // Compute the residuals
    Real maxRi = 0.0;

    for (int i = 0; i < lsp.m_numP*dim; i++)
    {
      Real AtimeX = 0.0;

      for (int j = 0; j < lsp.m_numP + lsp.m_numPLess1; j++)
      {
        AtimeX += lsp.m_matrix[i][j] * unknowns[j];
      }

      Real ri = AtimeX - rhs[i];

      if (Abs(ri) > maxRi)
      {
        m_residual[iDegree][0] = Abs(ri);
        maxRi = Abs(ri);
      }

      m_residual[iDegree][1] += Abs(ri);
      m_residual[iDegree][2] += ri * ri;
    }

    m_residual[iDegree][2] = sqrt(m_residual[iDegree][2]);
  }
}

template <int dim> inline void CutCellMoments<dim>::computeResiduals(const Vector< CutCellMoments<dim> > & a_refinedCCMoms,
                                                                     const int                           & a_degreeP)
{
  for (int iDegree = 0; iDegree <= a_degreeP; iDegree++)
  {
    Real maxRes = 0.0;
    Real tempRes2 = 0.0;
    Real tempRes1 = 0.0;

    for (int i = 0; i < a_refinedCCMoms.size(); i++)
    {
      Real res0 = a_refinedCCMoms[i].m_residual[iDegree][0];
      Real res1 = a_refinedCCMoms[i].m_residual[iDegree][1];
      Real res2 = a_refinedCCMoms[i].m_residual[iDegree][2];

      if (res0 > maxRes && res0 != LARGEREALVAL)
      {
        maxRes = res0;
        m_residual[iDegree][0] = res0;
      }

      if (res1 != LARGEREALVAL)
      {
        tempRes1 += res1;
      }

      if (res2 != LARGEREALVAL)
      {
        tempRes2 += res2*res2;
      }
    }

    m_residual[iDegree][1] = tempRes1;

    if (tempRes2 != LARGEREALVAL)
    {
      m_residual[iDegree][2] = sqrt(tempRes2);
    }
  }
}

// Methods for reading moments that check the moment exists
template <int dim> inline Real CutCellMoments<dim>::getMoment(const IvDim   & a_mono,
                                                              const EBorVol & a_EBorVol) const
{
  Real moment = LARGEREALVAL;

  if (a_EBorVol == VolMoment)
  {
    // Find the vol in the map
    typename PthMoment::const_iterator it = m_moments.find(a_mono);

    if (it != m_moments.end())
    {
      moment = it->second;
    }
    else
    {
      MayDay::Abort("No volume moment in m_moments");
    }
  }
  else if (a_EBorVol == EBMoment)
  {
    // Find the vol in the EBmap
    typename PthMoment::const_iterator it = m_EBmoments.find(a_mono);

    if (it != m_EBmoments.end())
    {
      moment = it->second;
    }
    else
    {
      MayDay::Abort("No volume moment in m_bdMoments");
    }
  }
  else
  {
    MayDay::Abort("Must ask for EBMoment or VolMoment");
  }

  return moment;
}

inline Real CutCellMoments<1>::getMoment(const IndexTM<int,1> & a_mono,
                                         const EBorVol        & a_EBorVOL) const
{
  return getMoment(a_mono);
}

// Method for reading moments
inline Real CutCellMoments<1>::getMoment(const IndexTM<int,1> & a_mono) const
{
  Real moment = LARGEREALVAL;

  // Find the vol in the map
  OneDMoments::const_iterator it = m_moments.find(a_mono);

  if (it != m_moments.end())
  {
    moment = it->second;
  }
  else
  {
    MayDay::Abort("No volume moment in top face moment map");
  }

  return moment;
}

template <int dim> inline Real CutCellMoments<dim>::getVol(const EBorVol& a_EBorVol) const
{
  Real volume;

  // Volume calculated by  integrating x^0 = 1
  IvDim zero;

  for (int idir = 0; idir < dim; ++idir)
  {
    zero[idir] = 0;
  }

  volume = getMoment(zero,a_EBorVol);

  return volume;
}

// Read geometric data and sanity check
inline Real CutCellMoments<1>::getVol(const EBorVol & a_EBorVol) const
{
  IndexTM<int,1> zero;
  zero[0] = 0;
  Real volume = getMoment(zero);

  return volume;
}

template <int dim> inline IndexTM<Real,dim> CutCellMoments<dim>::getCentroid(const EBorVol& a_EBorVol) const
{
  // This call checks the string to be either "VolMoment" or "EBMoment"
  Real volume = getVol(a_EBorVol);

  RvDim centroid;

  // Mono= 0,0,...1,0,0
  for (int idir = 0; idir < dim; ++idir)
  {
    IvDim mono     = BASISV_TM<int,dim>(idir);
    centroid[idir] = getMoment(mono,a_EBorVol);

    if (volume < 0.0)
    {
      MayDay::Warning("CutCellMoments::getCentroid: Volume fraction is negative");
    }
    else if (volume == 0.0)
    {
      MayDay::Error("CutCellMoments::getCentroid: Volume fraction is zero");
    }

    centroid[idir] /= volume;
  }

  return centroid;
}

inline IndexTM<Real,1> CutCellMoments<1>::getCentroid(const EBorVol& a_EBorVol) const
{
  IndexTM<Real,1> centroid;
  IndexTM<int,1> first;
  first[0] = 1;

  Real volume = getVol(a_EBorVol);

  if (volume <= 0)
  {
    MayDay::Abort("Volume equals 0");
  }
  else
  {
    centroid[0] = getMoment(first)/volume;
  }

  return centroid;
}

template <int dim> inline Real CutCellMoments<dim>::getResidual(const int & a_iDegree,
                                                                const int & a_normJ) const
{
  return m_residual[a_iDegree][a_normJ];
}

template <int dim> inline bool CutCellMoments<dim>::isCovered() const
{
  return m_iFData.m_allVerticesOut;
}

inline bool CutCellMoments<1>::isCovered() const
{
  return m_iFData.m_allVerticesOut;
}

template <int dim> inline bool CutCellMoments<dim>::isRegular() const
{
  return (m_iFData.m_allVerticesIn && (!m_bdCCOn));
}

inline bool CutCellMoments<1>::isRegular() const
{
  return m_iFData.m_allVerticesIn;
}

template <int dim> inline void CutCellMoments<dim>::print(ostream& a_out) const
{
   // Print m_pthMoment.
  a_out << "m_pthMoment:\n";

  for (typename PthMoment::const_iterator it = m_moments.begin(); it != m_moments.end(); ++it)
  {
    a_out << "  Integral " << it->first << " = " << it->second << '\n';
  }

  for (typename PthMoment::const_iterator it = m_EBmoments.begin(); it != m_EBmoments.end(); ++it)
  {
    a_out << "  EBIntegral " << it->first << " = " << it->second << '\n';
  }

  a_out << "IFData:\n" << m_iFData << "**" << "\n";

  a_out << "Dim = " << dim << '\n';
  a_out << "Number of bdCutCellMoments = " << m_bdCutCellMoments.size() << '\n';

  for (typename BdCutCellMoments::const_iterator it = m_bdCutCellMoments.begin(); it != m_bdCutCellMoments.end(); ++it)
  {
    a_out << "BdId = " << it->first << ", m_iFData = " << '\n';
    a_out << it->second.m_iFData << '\n';

    typedef map<IndexTM<int,dim-1>,Real,LexLT <IndexTM<int,dim-1> > > BdPthMoment;

    for (typename BdPthMoment::const_iterator itbd = it->second.m_moments.begin(); itbd != it->second.m_moments.end(); ++itbd)
    {
      a_out << "  Integral " << itbd->first << " = " << itbd->second << '\n';
    }

    for (typename BdPthMoment::const_iterator itbd = it->second.m_EBmoments.begin(); itbd != it->second.m_EBmoments.end(); ++itbd)
    {
      a_out << "  EBIntegral " << itbd->first << " = " << itbd->second << '\n';
    }
  }
}

inline void CutCellMoments<1>::print(ostream & a_out) const
{
  a_out << "IFData = " << m_iFData << "\n";
  for (OneDMoments::const_iterator it = m_moments.begin(); it != m_moments.end();++it)
  {
    a_out << "Moment: " << it->first << " = " << it->second << "\n";
  }
}

template <int dim> inline void CutCellMoments<dim>::dump() const
{
  print(pout());
}

template <int dim> inline void CutCellMoments<dim>::printMoments()
{
  pout() << "EB Moments" << "\n";
  for (typename PthMoment::const_iterator it = m_EBmoments.begin(); it != m_EBmoments.end(); ++it)
  {
    pout() << "EB Integral" << it->first << "=" << m_EBmoments[it->first] << "\n";
  }

  pout() << "Moments" << "\n";
  for (typename PthMoment::const_iterator it = m_moments.begin(); it != m_moments.end(); ++it)
  {
    pout() << "Integral" << it->first << "=" << m_moments[it->first] << "\n";
  }
}

template <int dim> inline void CutCellMoments<dim>::printBdMoments()
{
  IndexTM<int,2> bdId;

  for (int idir = 0; idir < dim; idir++)
  {
    bdId[0] = idir;

    for (int hilo = 0; hilo <= 1; hilo++)
    {
      bdId[1] = hilo;

      pout() << "bdId" << bdId << "\n";

      pout() << "EB Moments" << "\n";
      CutCellMoments<dim-1> bdCutCell = getBdCutCellMoments(bdId);
      for (typename PthMomentLesserDimension::const_iterator it = bdCutCell.m_EBmoments.begin(); it != bdCutCell.m_EBmoments.end(); ++it)
      {
        pout() << "EB Integral" << it->first << "=" << bdCutCell.m_EBmoments[it->first] << "\n";
      }

      pout() << "Moments" << "\n";
      for (typename PthMomentLesserDimension::const_iterator it = bdCutCell.m_moments.begin(); it != bdCutCell.m_moments.end(); ++it)
      {
        pout() << "Integral" << it->first << "=" << bdCutCell.m_moments[it->first] << "\n";
      }
    }
  }
}

// Operators
template <int dim> inline void CutCellMoments<dim>::operator=(const CutCellMoments<dim> & a_cutCellMoments)
{
  m_moments          = a_cutCellMoments.m_moments;
  m_EBmoments        = a_cutCellMoments.m_EBmoments;
  m_bdCutCellMoments = a_cutCellMoments.m_bdCutCellMoments;
  m_iFData           = a_cutCellMoments.m_iFData;
  m_bdCCOn           = a_cutCellMoments.m_bdCCOn;
  m_residual         = a_cutCellMoments.m_residual;
}

inline void CutCellMoments<1>::operator=(const CutCellMoments<1> & a_cutCellMoments)
{
  m_iFData  = a_cutCellMoments.m_iFData;
  m_moments = a_cutCellMoments.m_moments;
}

template <int dim> inline ostream& operator<<(ostream                   & a_out,
                                              const CutCellMoments<dim> & a_cutCellMoments)
{
  a_cutCellMoments.print(a_out);
  return a_out;
}

#include "NamespaceFooter.H"

#endif

---