MultiGrid.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: MultiGrid.hpp,v 1.13 2008/03/01 18:05:09 delpinux Exp $


#ifndef MULTIGRID_HPP
#define MULTIGRID_HPP

#include <Vector.hpp>

#include <Structured3DVector.hpp>
#include <Structured3DMeshShape.hpp>

#include <Structured3DMesh.hpp>

#include <SpectralMesh.hpp>

#include <UnAssembledMatrix.hpp>
#include <SparseMatrix.hpp>

#include <MultiGridOptions.hpp>
#include <GetParameter.hpp>

#include <Preconditioner.hpp>

#include <fstream>
#include <set>
#include <list>

#include <DegreeOfFreedomSet.hpp>

#include <BoundaryConditionDiscretizationElimination.hpp>

#include <SolverInformationCenter.hpp>

class MultiGrid
  : public Preconditioner
{
private:

  class DirichletPositions
  {
  private:
    Vector<size_t> __positions;

  public:
    size_t size() const
    {
      return __positions.size();
    }

    const size_t& operator()(const size_t i) const
    {
      ASSERT (i<__positions.size());
      return __positions[i];
    }

    void setPosition(const size_t& level, const size_t& pos)
    {
      __positions[level]=pos;
    }

    const DirichletPositions& operator=(const DirichletPositions& d)
    {
      __positions=d.__positions;
      return *this;
    }

    DirichletPositions(const size_t numberOfLevel)
      : __positions(numberOfLevel)
    {
      ;
    }

    DirichletPositions(const DirichletPositions& d)
      : __positions(d.__positions)
    {
      ;
    }

    ~DirichletPositions()
    {
      ;
    }
  };

  mutable size_t __level;
  int nu1;
  int nu2;
  int mu1;
  int mu2;

  int maxiter;

  real_t omega;
  real_t epsilon;

  const DegreeOfFreedomSet& __degreeOfFreedomSet;

  Structured3DMeshShape __meshShape;

  GetParameter<MultiGridOptions> __options;

  Structured3DVector<bool> __dirichlet;

  std::vector<Structured3DMeshShape> __meshesShape;
  std::list<DirichletPositions> __positionList;

  mutable std::vector<Vector<real_t> > __dirichletValues;

public:
  std::string name() const
  {
    return "finite differences multigrid";
  }


  MultiGrid(const Problem& problem,
            const SparseMatrix& A,
            const DegreeOfFreedomSet& degreeOfFreedomSet)
    : Preconditioner (problem,
                      Preconditioner::multigrid),
      __degreeOfFreedomSet(degreeOfFreedomSet),
      __meshShape(TinyVector<3,size_t>(0,0,0),
                  TinyVector<3,real_t>(-1,-1,-1),
                  TinyVector<3,real_t>(1,1,1)),
      __dirichlet(Array3DShape(0,0,0))
  {
    const Mesh& solverMesh = SolverInformationCenter::instance().mesh();
    const DiscretizationType& discretization = SolverInformationCenter::instance().discretizationType();
    switch (solverMesh.type()) {
    case Mesh::cartesianHexahedraMesh: {
      const Structured3DMesh& mesh = static_cast<const Structured3DMesh&>(solverMesh);
      __meshShape = Structured3DMeshShape(TinyVector<3,size_t>(mesh.shape().shape().nx(),
                                                               mesh.shape().shape().ny(),
                                                               mesh.shape().shape().nz()),
                                          mesh.shape().a(),mesh.shape().b());
      __dirichlet.resize(mesh.shape().shape());
      break;
    }
    case Mesh::spectralMesh: {
      const SpectralMesh& mesh = static_cast<const SpectralMesh&>(solverMesh);
      __meshShape = Structured3DMeshShape(TinyVector<3,size_t>(mesh.shape().shape().nx()-1,
                                                               mesh.shape().shape().ny()-1,
                                                               mesh.shape().shape().nz()-1),
                                          mesh.shape().a(),mesh.shape().b());
//       __meshShape = mesh.shape();
      __dirichlet.resize(__meshShape.shape());
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "Multigrid does not implemented for this kind of mesh",
                         ErrorHandler::normal);      
    }
    }

    __level = (int)(std::log(real_t(__meshShape.nx()-1))/std::log(2.));

    nu1   = __options.value().nu1();
    nu2   = __options.value().nu2();
    mu1   = __options.value().mu1();
    mu2   = __options.value().mu2();

    epsilon = __options.value().epsilon();
    maxiter = __options.value().maxiter();

    omega = __options.value().omega();

    Structured3DMesh mesh(__meshShape,
                          new VerticesCorrespondance(__meshShape.numberOfVertices()));

    ReferenceCounting<BoundaryConditionDiscretization>
      e = new BoundaryConditionDiscretizationElimination<Structured3DMesh,
                                                         ScalarDiscretizationTypeBase::lagrangianFEM1>(__problem,
                                                                                                       mesh,
                                                                                                       __degreeOfFreedomSet);

//     ReferenceCounting<BaseMatrix> a = new UnAssembledMatrix(0);
//     ReferenceCounting<BaseVector> b
//       = new Vector<real_t>(__degreeOfFreedomSet.size());

//     (*e).setSecondMember(a,b);

    for (size_t i=0; i<__degreeOfFreedomSet.size(); ++i) {
      __dirichlet[i] = (*e).dirichlet(i);
    }

//     __dirichlet = true;
//     for (size_t i=1; i<__dirichlet.nx()-1; ++i) {
//       for (size_t j=1; j<__dirichlet.ny()-1; ++j) {
//      for (size_t k=1; k<__dirichlet.nz()-1; ++k) {
//        __dirichlet(i,j,k) = false;

//      }
//       }
//     }    

    __meshesShape.reserve(__level);
    __dirichletValues.reserve(__level);

    for (size_t i=0; i<__level; ++i) {
      size_t nx = ((__meshShape.nx()-1)/(1<<i))+1;
      size_t ny = ((__meshShape.ny()-1)/(1<<i))+1;
      size_t nz = ((__meshShape.nz()-1)/(1<<i))+1;

      __meshesShape.push_back(Structured3DMeshShape(TinyVector<3,size_t>(nx,ny,nz),
                                                    __meshShape.a(),
                                                    __meshShape.b()));
    }

    for (size_t i=0; i<__level; ++i) {
      __dirichletValues
        .push_back(Vector<real_t>(__meshesShape[i].numberOfVertices()));
    }


    for (size_t i=0; i<__dirichlet.nx(); ++i) {
      for (size_t j=0; j<__dirichlet.ny(); ++j) {
        for (size_t k=0; k<__dirichlet.nz(); ++k) {
          if(__dirichlet(i,j,k)) {
            DirichletPositions p(__level);
            // Stores the first position.
            p.setPosition(0, __dirichlet.shape()(i,j,k));

            for(size_t l=1; l<__level; ++l) {
              const real_t x = i*__meshShape.hx();
              const real_t y = j*__meshShape.hy();
              const real_t z = k*__meshShape.hz();

              const size_t I = size_t(rint(x/__meshesShape[l].hx()));
              const size_t J = size_t(rint(y/__meshesShape[l].hy()));
              const size_t K = size_t(rint(z/__meshesShape[l].hz()));

              p.setPosition(l,__meshesShape[l].shape()(I,J,K));
            }
            __positionList.push_back(p);
          }
        }
      }
    }
  }

  void copyDirichlet(const Vector<real_t>& u,
                     Vector<real_t>& v,
                     const size_t level) const
  {
    ASSERT(v.size() == __meshesShape[__level-level].numberOfVertices());
    ASSERT(u.size() == __meshesShape[__level-level].numberOfVertices());

    for(std::list<DirichletPositions>::const_iterator i
          = __positionList.begin();
        i != __positionList.end(); ++i) {
      size_t j = (*i)(__level-level);
      v[j] = u[j];
    }
  }

  void residuDirichlet(Vector<real_t>& v,
                       const Vector<real_t>& u1,
                       const Vector<real_t>& u2,
                       const size_t level) const
  {
    ASSERT(v.size() == __meshesShape[__level-level].numberOfVertices());
    for(std::list<DirichletPositions>::const_iterator i
          = __positionList.begin();
        i != __positionList.end(); ++i) {
      size_t j = (*i)(__level-level);
      v[j] = u1[j] - u2[j];
    }
  }

  void setDirichlet(const size_t level,Vector<real_t>& v) const
  {
    for(std::list<DirichletPositions>::const_iterator i
          = __positionList.begin();
        i != __positionList.end(); ++i) {
      size_t j = (*i)(__level-level);
      v[j] = __dirichletValues[__level-level][j];
    }
  }

  void getDirichlet(const size_t level, Vector<real_t>& v) const
  {
    for(std::list<DirichletPositions>::const_iterator i
          = __positionList.begin();
        i != __positionList.end(); ++i) {
      size_t j = (*i)(__level-level);
      __dirichletValues[__level-level][j] = v[j];
    }
  }

  void interpolateDirichlet(Vector<real_t>& u2,
                            const Vector<real_t>& u1,
                            const size_t currentLevel) const
  {
    for(std::list<DirichletPositions>::const_iterator i
          = __positionList.begin();
        i != __positionList.end(); ++i) {
      
      const size_t i1 = (*i)(__level-currentLevel);
      const size_t i2 = (*i)(__level-currentLevel+1);

      u2[i2] = u1[i1];
    }
    //    getDirichlet(currentLevel-1,u1);
  }

  void restrictDirichlet(const Vector<real_t>& u2,
                         Vector<real_t>& u1,
                         const size_t currentLevel) const
  {
    for(std::list<DirichletPositions>::const_iterator i
          = __positionList.begin();
        i != __positionList.end(); ++i) {
      
      const size_t i1 = (*i)(__level-currentLevel);
      const size_t i2 = (*i)(__level-currentLevel+1);

      u1[i1] = u2[i2];
    }
  }

  void initializes()
  {
    ;
  }

  void weightJacobi(Structured3DVector<real_t>& u,
                    const Structured3DVector<real_t>& b,
                    const size_t currentLevel) const
  {
    const Array3DShape& shape = b.shape();

    Vector<real_t> d (u.size());

    real_t dx= __meshShape.dimension(0)/(shape.nx()-1);
    real_t dy= __meshShape.dimension(1)/(shape.ny()-1);
    real_t dz= __meshShape.dimension(2)/(shape.nz()-1);

    Structured3DVector<real_t> v(u.shape());
    v=0;

    for (size_t i=1; i<shape.nx()-1; ++i)
      for (size_t j=1; j<shape.ny()-1; ++j)
        for (size_t k=1; k<shape.nz()-1; ++k) {
          real_t dijk = (2./(dx*dx)+2./(dy*dy)+2./(dz*dz)
//                         +1E5*(pow((i*dx-0.5),2)+pow((j*dy-0.5),2)+pow((k*dz-0.5),2)<0.0625)
                         );
          v(i,j,k) = ((u(i-1,j,k) + u(i+1,j,k))/(dx*dx)
                      + (u(i,j-1,k) + u(i,j+1,k))/(dy*dy)
                      + (u(i,j,k-1) + u(i,j,k+1))/(dz*dz)
                      + b(i,j,k)
                      )/dijk;
        }

    copyDirichlet(u,v,currentLevel);

    v *= 1-omega;
    u *= omega;
    u += v;
  }

  void computeResidu(Structured3DVector<real_t>& r,
                     const Structured3DVector<real_t>& u,
                     const Structured3DVector<real_t>& f,
                     const size_t currentLevel) const
  {
    real_t dx= __meshShape.dimension(0)/(u.nx()-1);
    real_t dy= __meshShape.dimension(1)/(u.ny()-1);
    real_t dz= __meshShape.dimension(2)/(u.nz()-1);

    r = 0.;

    ASSERT (u.shape()==r.shape());
    ASSERT (u.shape()==f.shape());

      for (size_t i=1; i<u.nx()-1; ++i)
        for (size_t j=1; j<u.ny()-1; ++j)
          for (size_t k=1; k<u.nz()-1; ++k)
            r(i,j,k) = 
              f(i,j,k)
              - (u(i,j,k)*(2./(dx*dx)+2./(dy*dy)+2./(dz*dz)
//                           +1E5*(pow((i*dx-0.5),2)
//                                 +pow((j*dy-0.5),2)
//                                 +pow((k*dz-0.5),2)<0.0625)
                           )
                 - (u(i-1,j,k) + u(i+1,j,k))/(dx*dx)
                 - (u(i,j-1,k) + u(i,j+1,k))/(dy*dy)
                 - (u(i,j,k-1) + u(i,j,k+1))/(dz*dz));

      residuDirichlet(r,f,u,currentLevel);
  }


  void correct(Structured3DVector<real_t>& u,
               const Structured3DVector<real_t>& u2,
               const size_t currentLevel) const
  {

    getDirichlet(currentLevel,u);
    ASSERT(u.nx() == ((u2.nx()-1)*2+1));
    ASSERT(u.ny() == ((u2.ny()-1)*2+1));
    ASSERT(u.nz() == ((u2.nz()-1)*2+1));

    for (size_t i=1; i<u.nx()-1; ++i) {
      std::vector<int> indicesI;
      indicesI.push_back(i/2);
      if (i%2 != 0)
        indicesI.push_back(i/2 + 1);
      for (size_t j=1; j<u.ny()-1; ++j) {
        std::vector<int> indicesJ;
        indicesJ.push_back(j/2);
        if (j%2 != 0)
          indicesJ.push_back(j/2 + 1);
        for (size_t k=1; k<u.nz()-1; ++k) {
          std::vector<int> indicesK;
          indicesK.push_back(k/2);
          if (k%2 != 0)
            indicesK.push_back(k/2 + 1);

          real_t sum=0;
          for (size_t ii=0; ii<indicesI.size(); ++ii) {
            for (size_t jj=0; jj<indicesJ.size(); ++jj) {
              for (size_t kk=0; kk<indicesK.size(); ++kk) {
                sum += u2(indicesI[ii],indicesJ[jj],indicesK[kk]);
              }
            }
          }
          u(i,j,k)
            += sum / (indicesI.size()*indicesJ.size()*indicesK.size());
        }
      }
    }
    setDirichlet(currentLevel,u);
    //    restrictDirichlet(u2,u,currentLevel);
  }

  void project (Structured3DVector<real_t>& u2,
                const Structured3DVector<real_t>& u,
                const size_t currentLevel) const
  {
    //    u2=0;

    real_t c0 = 1./8.;
    real_t c1 = 1./16;
    real_t c2 = 1./32;
    real_t c3 = 1./64;

    for (size_t i=1; i<u2.nx()-1; ++i)
      for (size_t j=1; j<u2.ny()-1; ++j)
        for (size_t k=1; k<u2.nz()-1; ++k) {
          int i1 = 2*i;
          int i2 = 2*j;
          int i3 = 2*k;

          u2(i,j,k) = c0 * u(i1,i2,i3)
            + c1 * ( u(i1-1,i2,  i3  ) + u(i1+1,i2,  i3  )
                     +u(i1,  i2-1,i3  ) + u(i1,  i2+1,i3  )
                     +u(i1,  i2,  i3-1) + u(i1,  i2,  i3+1) )
            + c2 * ( u(i1-1,i2-1,i3  ) + u(i1+1,i2-1,i3  )
                     +u(i1-1,i2+1,i3  ) + u(i1+1,i2+1,i3  )
                     +u(i1,  i2-1,i3-1) + u(i1,  i2+1,i3-1)
                     +u(i1,  i2-1,i3+1) + u(i1,  i2+1,i3+1)
                     +u(i1-1,i2,  i3-1) + u(i1-1,i2,  i3+1)
                     +u(i1+1,i2,  i3-1) + u(i1+1,i2,  i3+1) )
            + c3 * ( u(i1-1,i2-1,i3-1) + u(i1+1,i2-1,i3-1)
                     +u(i1-1,i2+1,i3-1) + u(i1+1,i2+1,i3-1)
                     +u(i1-1,i2-1,i3+1) + u(i1+1,i2-1,i3+1)
                     +u(i1-1,i2+1,i3+1) + u(i1+1,i2+1,i3+1) );
        }
    interpolateDirichlet(u2,u,currentLevel);
  }

  void multigrid(Structured3DVector<real_t>& u,
                 const Structured3DVector<real_t>& b,
                 int currentLevel) const
  {
    Structured3DVector<real_t> r(u.shape());

    for (int i=0; i<nu1; ++i)
      weightJacobi(u,b,currentLevel);

    Array3DShape newShape ((u.nx()+1)/2,
                           (u.ny()+1)/2,
                           (u.nz()+1)/2);

    if (currentLevel > 1) {
      
      for (int i=0; i<mu1; ++i) { // mu-Cycle
        Structured3DVector<real_t> u2(newShape);
        u2 = 0;

        Structured3DVector<real_t> b2(newShape);
        computeResidu(r,u,b, currentLevel);
        
        project(b2,r,currentLevel);

        multigrid(u2, b2, currentLevel-1);
        correct(u,u2,currentLevel);
      }
    } else {
#warning Arrange this.
      for (size_t i=0; i<nu1+nu2; ++i)
        weightJacobi(u,b,currentLevel);
    }
    for (int i=0; i<nu2; ++i)
      weightJacobi(u,b,currentLevel);
  }

  void FMV(Structured3DVector<real_t>& u,
           const Structured3DVector<real_t>& b,
           int currentLevel) const
  {
    if (currentLevel > 1) {
      Array3DShape newShape ((u.nx()+1)/2,
                             (u.ny()+1)/2,
                             (u.nz()+1)/2);

      Structured3DVector<real_t> r(u.shape());

      Structured3DVector<real_t> u2(newShape);
      u2 = 0;

      Structured3DVector<real_t> b2(newShape);
      computeResidu(r,u,b,currentLevel);
      project(b2,r,currentLevel);

      FMV(u2, b2, currentLevel-1);
      correct(u,u2,currentLevel);
    }

    for (int i=0; i<mu2; ++i)
      multigrid(u,b,currentLevel);
  }

  void computes(const Vector<real_t>& f,
                Vector<real_t>& v) const
  {
    Structured3DVector<real_t> u(__meshShape.shape());
    Structured3DVector<real_t> b(__meshShape.shape());

    ASSERT(u.size() == v.size());
    ASSERT(b.size() == f.size());

    u = v;
    b = f;

    real_t dx = __meshShape.hx();
    real_t dy = __meshShape.hy();
    real_t dz = __meshShape.hz();

//     std::cout << "multigrid dxi = "<< dx << 'x' << dy << 'x' << dz << '\n';

//     real_t pi = std::atan(1.)*4.;

//     for (size_t i=0; i<b.nx(); ++i)
//       for (size_t j=0; j<b.ny(); ++j)
//      for (size_t k=0; k<b.nz(); ++k) {
//        b(i,j,k) = std::sin(pi*i*__meshShape.hx());
//      }
//     v = b;
//     for (size_t i=1; i<b.nx()-1; ++i)
//       for (size_t j=1; j<b.ny()-1; ++j)
//      for (size_t k=1; k<b.nz()-1; ++k) {
//        real_t x = i*__meshShape.hx();
//        b(i,j,k) = pi*pi*std::sin(pi*x);
//      }

    getDirichlet(__level,b);

     //     b *= 1./(dx*dy*dz);

    // inside the domain 2nd member should be rescaled
    // for compatibility between FEM and FDM.
//     for (size_t i=1; i<b.nx()-1; ++i)
//       for (size_t j=1; j<b.ny()-1; ++j)
//      for (size_t k=1; k<b.nz()-1; ++k) {
//        b(i,j,k) /= dx*dy*dz;
//        //      u(i,j,k) = 0;
//      }

//     for (size_t i=1; i<b.nx()-1; ++i)
//       for (size_t j=1; j<b.ny()-1; ++j)
//      for (size_t k=1; k<b.nz()-1; ++k) {
//        u(i,j,k) = 0;
//      }

//     for (size_t i=1; i<b.nx()-1; ++i)
//       for (size_t j=1; j<b.ny()-1; ++j)
//      for (size_t k=1; k<b.nz()-1; ++k) {
//        b(i,j,k) /= dx*dy*dz*dx*dy*dz;
//      }
//    b *= (dx*dy*dz)*(dx*dy*dz);
    u=0;

    setDirichlet(__level,u);
    setDirichlet(__level,b);

    Structured3DVector<real_t> r(u);
    computeResidu(r,u,b,__level);

    real_t resid0 = Norm(r);

    real_t norm;
    int n = 0;
    /*
    for (size_t i=0; i<10; ++i) {
      weightJacobi(u,b,__level);

      computeResidu(r,u,b,__level);
      real_t resid = Norm(r);

      ffout(4) << "\tFMG norme:" << resid << '\n';
    }
    */
//     computeResidu(r,u,b,__level);
//     ffout(0) << "\tFMG norme:" << Norm(r) << "\n";

    do {
      FMV(u, b, __level);

      computeResidu(r,u,b,__level);

      real_t resid = Norm(r);

      norm = resid/resid0;
      ffout(0) << "\tFMG norme:" << norm <<" (" << resid << '/' << resid0 << ")\n";
      n++;
    } while ((norm>epsilon) and (n<maxiter));

    for (size_t i=1; i<b.nx()-1; ++i)
      for (size_t j=1; j<b.ny()-1; ++j)
        for (size_t k=1; k<b.nz()-1; ++k) {
          v[u.shape()(i,j,k)] = u(i,j,k);
//        b(i,j,k) /= dx*dy*dz*dx*dy*dz;
        }
//     std::cout << "Norm(u)=" << Norm(u) << '\n';
  }

//   void Iterate(const Vector<real_t>& b,
//             SparseMatrix<MatType>&A,
//             Vector<real_t>& u)
//   {
//     Vector<MatType> D(A.size());
//     for(size_t i=0; i<A.size();++i)
//       D(i) = A(i,i);
//     SparseMatrix<MatType> E(A.size(), 13);
//     SparseMatrix<MatType> F(A.size(), 13);
//     for (size_t i=0; i<A.size(); ++i) {
//       const SparseLine<MatType>& Ai = A.Line(i);
//       for (size_t j=0; j<Ai.size(); ++j) {
//      if (Ai[j].Column() == -1)
//        break;
//      if (Ai[j].Column() < i)
//        E(i,Ai[j].Column()) = -Ai[j].Value();
//      if (Ai[j].Column() > i)
//        F(i,Ai[j].Column()) = -Ai[j].Value();
//       }
//     }

//     Vector<real_t> u1(u.size());
//     u1 = E*u + F*u;

//     for (size_t i=0; i<u1.size(); ++i)
//       u1[i] /= D[i];

//     u1 *= omega;
//     u  *= (1-omega);

//     u += u1;

//     u1 = b;
//     for (size_t i=0; i<u1.size(); ++i)
//       u1[i] /= D[i];
//     u1 *= omega;

//     u += u1;


//     //! Computes the multigrid interpolation.
    

//     vector<int> newShape(3);
//     newShape[0] = (shape[0]+1)/2;
//     newShape[1] = (shape[1]+1)/2;
//     newShape[2] = (shape[2]+1)/2;

//     SparseMatrix<real_t> IhH(u2.size(), 27);
//     for (size_t i=0; i<IhH.size(); ++i)
//       IhH(i,2*i) = 1.; // Dirichlet boundary conditions

//     for (size_t i=1; i<newShape[0]-1; ++i) {
//       for (size_t j=1; j<newShape[1]-1; ++j) {
//      for (size_t k=1; k<newShape[2]-1; ++k) {
//        int ijk   = i*newShape[1]*newShape[2]+j*newShape[2]+k;
//        int IJK   = (i*2)*shape[1]*shape[2]+(j*2)*shape[2]+2*k;
//        int I_1JK = (i*2-1)*shape[1]*shape[2]+(j*2)*shape[2]+2*k;
//        int I1JK  = (i*2+1)*shape[1]*shape[2]+(j*2)*shape[2]+2*k;
//        int IJ_1K = (i*2-1)*shape[1]*shape[2]+(j*2-1)*shape[2]+2*k;
//        int IJ1K  = (i*2+1)*shape[1]*shape[2]+(j*2+1)*shape[2]+2*k;
//        int IJK_1 = (i*2-1)*shape[1]*shape[2]+(j*2)*shape[2]+2*k-1;
//        int IJK1  = (i*2+1)*shape[1]*shape[2]+(j*2)*shape[2]+2*k+1;
//        IhH(ijk, IJK) = 1./2;
//        IhH(ijk, I_1JK) = 1./12;
//        IhH(ijk, I1JK) = 1./12;
//        IhH(ijk, IJ_1K) = 1./12;
//        IhH(ijk, IJ1K) = 1./12;
//        IhH(ijk, IJK_1) = 1./12;
//        IhH(ijk, IJK1) = 1./12;
//      }
//       }
//     }

//     // Coarse grid vectors ...
//     Vector<real_t> uH(newShape[0]*newShape[1]*newShape[2]);
//     Vector<real_t> bH(newShape[0]*newShape[1]*newShape[2]);

//     // Coarse matrix
//     SparseMatrix<MatType> AH(uH.size(), 27);

//     uH = IhH*u;
//     bH = IhH*b;
//     for (size_t i=0; i<u2.size(); ++i)
//       fferr(2) << u2[i] << '\n';

//   }

};
#endif // _MULTIGRID_HPP_

---