SpectralFEMPreconditioner.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stephane 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: SpectralFEMPreconditioner.cpp,v 1.19 2008/03/10 21:58:10 delpinux Exp $

#include <SpectralFEMPreconditioner.hpp>

#include <MultiGrid.hpp>

#include <Structured3DMesh.hpp>
#include <SpectralMesh.hpp>

#include <DiscretizationType.hpp>
#include <DegreeOfFreedomSetBuilder.hpp>

#include <MatrixManagement.hpp>
#include <SolverInformationCenter.hpp>

#include <FEMDiscretization.hpp>
#include <BoundaryConditionDiscretizationFEM.hpp>

#include <ScalarDiscretizationTypeSpectral.hpp>

#include <SpectralConformTransformation.hpp>

#include <LegendreBasis.hpp>
#include <LagrangeBasis.hpp>

#include <IdentityPrecond.hpp>
#include <DiagPrecond.hpp>
#include <ConjugateGradient.hpp>

#include <FEMFunctionBuilder.hpp>
#include <FEMFunction.hpp>

#include <VariationalProblem.hpp>
#include <ScalarFunctionConstant.hpp>
#include <VariationalOperatorAlphaUV.hpp>
#include <VariationalOperatorFV.hpp>

#include <SpectralLegendreDiscretizationConform.hpp>
#include <BoundaryConditionDiscretizationSpectralConform.hpp>

#include <SparseMatrix.hpp>

class SpectralFEMPreconditioner::Internal
{
private:
  const Problem& __problem;

  mutable TinyVector<3,size_t> __legendreDegrees;
  mutable TinyVector<3,size_t> __legendreDimensions;

  mutable TinyVector<3,size_t> __lagrangeDegrees;
  mutable TinyVector<3,size_t> __lagrangeDimensions;

  mutable TinyVector<3,real_t> __a;
  mutable TinyVector<3,real_t> __b;

  mutable TinyVector<3,Vector<real_t> > __lagrangeVertices;

  template <typename MeshType>
  void __solveFEM(ConstReferenceCounting<MeshType> mesh,
                  Structured3DVector<real_t>& lagrange) const
  {
    ParameterCenter::instance().set("memory::matrix", "sparse");

    if (__problem.type() != Problem::variational) {
      throw ErrorHandler(__FILE__,__LINE__,
                         "not yet implemented",
                         ErrorHandler::unexpected);
    }

    VariationalProblem femProblem(dynamic_cast<const VariationalProblem&>(__problem), true);

    FEMFunctionBuilder builder;
    builder.build(ScalarDiscretizationTypeFEM(ScalarDiscretizationTypeBase::lagrangianFEM1),
                  mesh,
                  lagrange);

    femProblem.add(new VariationalOperatorFV(0,VariationalOperator::normal, builder.getBuiltScalarFunction()));
    
    DiscretizationType discretization;
    for (size_t i=0; i<femProblem.numberOfUnknown(); ++i) {
      discretization.add(new ScalarDiscretizationTypeFEM(ScalarDiscretizationTypeBase::lagrangianFEM1));
    }

    DegreeOfFreedomSetBuilder dofBuilder(discretization,
                                         *mesh);
    const DegreeOfFreedomSet& degreeOfFreedomSet
      = dofBuilder.degreeOfFreedomSet();

    MemoryManager MM;
    ReferenceCounting<BaseMatrix> femA;

    bool performAssembling =MM.ReserveMatrix(femA,
                                             femProblem.numberOfUnknown(),
                                             degreeOfFreedomSet.size());

//     ReferenceCounting<BaseVector> femB;
//     MM.ReserveVector(femB,
//                   femProblem.numberOfUnknown(),
//                   degreeOfFreedomSet.size());

//     ffout(0) << "Finite element method: disretization...\n";

//     ReferenceCounting<FEMDiscretization<MeshType,
//       ScalarDiscretizationTypeBase::lagrangianFEM1> > FEM
//       = new FEMDiscretization<MeshType,
//       ScalarDiscretizationTypeBase::lagrangianFEM1>(femProblem,
//                                                  *mesh,
//                                                  discretization,
//                                                  *femA, *femB, degreeOfFreedomSet);

//     if (performAssembling) {
//       FEM->assembleMatrix();
//     } else {
//       ffout(4) << "- keeping previous operator discretization\n";
//     }

//     FEM->assembleSecondMember();


//     ffout(4) << "- discretizing boundary conditions\n";

//     BoundaryConditionDiscretizationFEM<MeshType,
//       ScalarDiscretizationTypeBase::lagrangianFEM1>* bcd
//       = new BoundaryConditionDiscretizationFEM<MeshType,
//       ScalarDiscretizationTypeBase::lagrangianFEM1>(femProblem,
//                                                  *mesh,
//                                                  degreeOfFreedomSet);
//     bcd->associatesMeshesToBoundaryConditions();
//     ReferenceCounting<BoundaryConditionDiscretization> bcDiscretization = bcd;

//     // Set Dirichlet information to the matrix
//     FEM->setDirichletList(bcDiscretization->getDirichletList());

//     ffout(4) << "- second member modification\n";
//     bcDiscretization->setSecondMember(femA,femB);

//     ffout(4) << "- matrix modification\n";
//     bcDiscretization->setMatrix(femA,femB);

//     ffout(4) << "Finite element method: disretization done\n";

    if (femA->type() == BaseMatrix::doubleHashedMatrix) {
      Timer t;
      t.start();

      SparseMatrix* aa
        = new SparseMatrix(static_cast<DoubleHashedMatrix&>(*femA));
    
      femA = aa; // now use sparse matrix

      t.stop();
      ffout(4) << "Matrix copy: " << t << '\n';
    }

//     ReferenceCounting<Preconditioner> P = new IdentityPrecond(femProblem);

//     lagrange = 0;

    Structured3DVector<real_t> residu(lagrange);
    MultiGrid M(__problem,static_cast<const SparseMatrix&>(*femA),degreeOfFreedomSet);
    M.computes(residu,lagrange);
    //    static_cast<Vector<real_t>&>(*femB)=0;
//     ffout(4) << "- FEM CG\n";
//     ConjugateGradient cg(static_cast<Vector<real_t>&>(*femB),
//                       *femA, *P,
//                       static_cast<Vector<real_t>&>(lagrange),500,1E-2);
//     //    lagrange = static_cast<Vector<real_t>&>(*femB);
//     ffout(4) << "- FEM CG: done\n"; 
  }


  void __legendreToLagrange(const Structured3DVector<real_t>& legendre,
                            Structured3DVector<real_t>& lagrange) const
  {
    lagrange = 0;
    
    TinyVector<3, ConstReferenceCounting<LegendreBasis> > legendreBasis;
    for (size_t i=0; i<3; ++i) {
      legendreBasis[i] = new LegendreBasis(__legendreDegrees[i]);
    }

    TinyVector<3,Vector<Vector<real_t> > >  legendreBaseValues;
    for (size_t i=0; i<3; ++i) {
      legendreBaseValues[i].resize(__lagrangeVertices[i].size());
    }

    for(size_t i=0; i<3; ++i) {
      const Vector<real_t> vertices = __lagrangeVertices[i];
      for(size_t j=0; j<vertices.size(); ++j){
        // Compute position of vertex j in direction i
        const real_t& xj = vertices[j];

        legendreBaseValues[i][j].resize(legendreBasis[i]->dimension());
        legendreBasis[i]->getValues(xj,legendreBaseValues[i][j]);
      }
    }

    const Structured3DVector<real_t>& u_ijk = legendre;

    Structured3DVector<real_t> u_ijt(Array3DShape(__legendreDimensions[0],
                                                  __legendreDimensions[1],
                                                  __lagrangeDimensions[2]));
    u_ijt=0;

    {
      for (size_t i=0; i< __legendreDimensions[0]; ++i) {
        for (size_t j=0; j< __legendreDimensions[1]; ++j) {
          for (size_t k=0; k< __legendreDimensions[2]; ++k) {
            const real_t& u_ijk_value = u_ijk(i,j,k);
            for (size_t t=0; t< __lagrangeDimensions[2]; ++t) {
              u_ijt(i,j,t) += u_ijk_value * legendreBaseValues[2][t][k];
            }
          }
        }
      }
    }

    Structured3DVector<real_t> u_ist(Array3DShape(__legendreDimensions[0],
                                                  __lagrangeDimensions[1],
                                                  __lagrangeDimensions[2]));
    u_ist=0;

    for (size_t i=0; i<__legendreDimensions[0]; ++i) {
      for (size_t j=0; j<__legendreDimensions[1]; ++j) {
        for (size_t t=0; t<__lagrangeDimensions[2]; ++t) {
          const real_t& u_ijt_value = u_ijt(i,j,t);
          for (size_t s=0; s<__lagrangeDimensions[1]; ++s) {
            u_ist(i,s,t) += u_ijt_value*legendreBaseValues[1][s][j];
          }
        }
      }
    }

    Structured3DVector<real_t>& u_rst = lagrange;
    u_rst=0;
    
    for (size_t i=0; i<__legendreDimensions[0]; ++i) {
      for (size_t s=0; s<__lagrangeDimensions[1]; ++s) {
        for (size_t t=0; t<__lagrangeDimensions[2]; ++t) {
          const real_t& u_ist_value = u_ist(i,s,t);
          for (size_t r=0; r<__lagrangeDimensions[0]; ++r) {
            u_rst(r,s,t) += u_ist_value*legendreBaseValues[0][r][i];
          }
        }
      }
    }
  }
  
  void __legendreToLagrangeTransposed(const Structured3DVector<real_t>& lagrange,
                                      Structured3DVector<real_t>& legendre) const
  {
    legendre = 0;
    
    TinyVector<3, ConstReferenceCounting<LegendreBasis> > legendreBasis;
    for (size_t i=0; i<3; ++i) {
      legendreBasis[i] = new LegendreBasis(__legendreDegrees[i]);
    }

    TinyVector<3,Vector<Vector<real_t> > >  legendreBaseValues;
    for (size_t i=0; i<3; ++i) {
      legendreBaseValues[i].resize(__lagrangeVertices[i].size());
    }

    for(size_t i=0; i<3; ++i) {
      const Vector<real_t> vertices = __lagrangeVertices[i];
      for(size_t j=0; j<vertices.size(); ++j){
        // Compute position of vertex j in direction i
        const real_t& xj = vertices[j];

        legendreBaseValues[i][j].resize(legendreBasis[i]->dimension());
        legendreBasis[i]->getValues(xj,legendreBaseValues[i][j]);
      }
    }

    const Structured3DVector<real_t>& u_rst = lagrange;

    Structured3DVector<real_t> u_rsk(Array3DShape(__lagrangeDimensions[0],
                                                  __lagrangeDimensions[1],
                                                  __legendreDimensions[2]));
    u_rsk=0;

    {
      for (size_t r=0; r< __lagrangeDimensions[0]; ++r) {
        for (size_t s=0; s< __lagrangeDimensions[1]; ++s) {
          for (size_t t=0; t< __lagrangeDimensions[2]; ++t) {
            const real_t& u_rst_value = u_rst(r,s,t);
            for (size_t k=0; k< __legendreDimensions[2]; ++k) {
              u_rsk(r,s,k) += u_rst_value * legendreBaseValues[2][t][k];
            }
          }
        }
      }
    }

    Structured3DVector<real_t> u_rjk(Array3DShape(__lagrangeDimensions[0],
                                                  __legendreDimensions[1],
                                                  __legendreDimensions[2]));
    u_rjk=0;

    for (size_t r=0; r<__lagrangeDimensions[0]; ++r) {
      for (size_t s=0; s<__lagrangeDimensions[1]; ++s) {
        for (size_t k=0; k<__legendreDimensions[2]; ++k) {
          const real_t& u_rsk_value = u_rsk(r,s,k);
          for (size_t j=0; j<__legendreDimensions[1]; ++j) {
            u_rjk(r,j,k) += u_rsk_value*legendreBaseValues[1][s][j];
          }
        }
      }
    }

    Structured3DVector<real_t>& u_ijk = legendre;
    u_ijk=0;
    
    for (size_t r=0; r<__lagrangeDimensions[0]; ++r) {
      for (size_t j=0; j<__legendreDimensions[1]; ++j) {
        for (size_t k=0; k<__legendreDimensions[2]; ++k) {
          const real_t& u_rjk_value = u_rjk(r,j,k);
          for (size_t i=0; i<__legendreDimensions[0]; ++i) {
            u_ijk(i,j,k) += u_rjk_value*legendreBaseValues[0][r][i];
          }
        }
      }
    }
  }

  void __lagrangeToLegendre(const Structured3DVector<real_t>& lagrange,
                            Structured3DVector<real_t>& legendre) const
  {
    legendre = 0;
    
    TinyVector<3, ConstReferenceCounting<LegendreBasis> > legendreBasis;
    TinyVector<3, ConstReferenceCounting<LagrangeBasis> > lagrangeBasis;
    for (size_t i=0; i<3; ++i) {
      legendreBasis[i] = new LegendreBasis(__legendreDegrees[i]);
      lagrangeBasis[i] = new LagrangeBasis(__lagrangeVertices[i]);
    }
    
    TinyVector<3, ConstReferenceCounting<GaussLobatto> > gaussLobatto;
    for (size_t i=0; i<3; ++i) {
      gaussLobatto[i] = GaussLobattoManager::instance().getReference((__legendreDegrees[i]+__lagrangeDegrees[i]) + 1);
    }

    TinyVector<3,Vector<Vector<real_t> > >  legendreBaseValues;
    TinyVector<3,Vector<Vector<real_t> > >  lagrangeBaseValues;
    for (size_t i=0; i<3; ++i) {
      legendreBaseValues[i].resize(gaussLobatto[i]->numberOfPoints());
      lagrangeBaseValues[i].resize(gaussLobatto[i]->numberOfPoints());
    }
    
    for(size_t i=0; i<3; ++i) {
      const GaussLobatto& gaussLobattoPoints = *gaussLobatto[i];
      for(size_t j=0; j<gaussLobattoPoints.numberOfPoints(); ++j){
        
        // Compute position of vertex j in direction i
        const real_t& xj = gaussLobattoPoints(j);
        
        legendreBaseValues[i][j].resize(legendreBasis[i]->dimension());
        legendreBasis[i]->getValues(xj,legendreBaseValues[i][j]);
        
        lagrangeBaseValues[i][j].resize(lagrangeBasis[i]->dimension());
        lagrangeBasis[i]->getValues(xj,lagrangeBaseValues[i][j]);
      }
    }
    
    Structured3DVector<real_t> u_r0s0t(Array3DShape(__lagrangeDimensions[0],
                                                    __lagrangeDimensions[1],
                                                    gaussLobatto[2]->numberOfPoints()));
    u_r0s0t = 0;
    
    {
      for (size_t r0=0; r0 < __lagrangeDimensions[0]; ++r0) {
        for (size_t s0=0; s0 < __lagrangeDimensions[1]; ++s0) {
          for (size_t t0=0; t0 < __lagrangeDimensions[2]; ++t0) {
            const real_t& u_r0s0t0_value = lagrange(r0,s0,t0);
            for (size_t t=0; t< gaussLobatto[2]->numberOfPoints(); ++t) {
              u_r0s0t(r0,s0,t) += u_r0s0t0_value * lagrangeBaseValues[2][t][t0] ;
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_r0st(Array3DShape(__lagrangeDimensions[0],
                                                   gaussLobatto[1]->numberOfPoints(),
                                                   gaussLobatto[2]->numberOfPoints()));
    u_r0st = 0;
    
    {
      for (size_t r0=0; r0 <__lagrangeDimensions[0]; ++r0) {
        for (size_t s0=0; s0 <__lagrangeDimensions[1]; ++s0) {
          for (size_t t=0; t < gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& u_r0s0t_value = u_r0s0t(r0,s0,t) ;
            for (size_t s=0; s<gaussLobatto[1]->numberOfPoints(); ++s) {
              u_r0st(r0,s,t) += u_r0s0t_value * lagrangeBaseValues[1][s][s0];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_rst(Array3DShape(gaussLobatto[0]->numberOfPoints(),
                                                  gaussLobatto[1]->numberOfPoints(),
                                                  gaussLobatto[2]->numberOfPoints()));
    u_rst = 0;
    {
      for (size_t r0=0; r0 < __lagrangeDimensions[0]; ++r0) {
        for (size_t s=0; s <gaussLobatto[1]->numberOfPoints(); ++s) {
          for (size_t t=0; t <gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& u_r0st_value = u_r0st(r0,s,t);
            for (size_t r=0; r<gaussLobatto[0]->numberOfPoints(); ++r) {
              u_rst(r,s,t) += u_r0st_value * lagrangeBaseValues[0][r][r0];
            }
          }
        }
      }
    }

    Structured3DVector<real_t> u_rsk(Array3DShape(gaussLobatto[0]->numberOfPoints(),
                                                  gaussLobatto[1]->numberOfPoints(),
                                                  __legendreDimensions[2]));
    u_rsk = 0;
    
    {
      for (size_t r=0; r < gaussLobatto[0]->numberOfPoints(); ++r) {
        for (size_t s=0; s < gaussLobatto[1]->numberOfPoints(); ++s) {
          for (size_t t=0; t < gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& wt = gaussLobatto[2]->weight(t);
            const real_t u_rst_value = u_rst(r,s,t) * wt;
            for (size_t k=0; k<__legendreDimensions[2]; ++k) {
              u_rsk(r,s,k) += u_rst_value * legendreBaseValues[2][t][k];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_rjk(Array3DShape(gaussLobatto[0]->numberOfPoints(),
                                                  __legendreDimensions[1],
                                                  __legendreDimensions[2]));
    u_rjk = 0;
    
    {
      for (size_t r=0; r < gaussLobatto[0]->numberOfPoints(); ++r) {
        for (size_t s=0; s < gaussLobatto[1]->numberOfPoints(); ++s) {
          const real_t& ws = gaussLobatto[1]->weight(s);
          for (size_t k=0; k < __legendreDimensions[2]; ++k) {
            const real_t u_rsk_value = u_rsk(r,s,k) * ws;
            for (size_t j=0; j<__legendreDimensions[1]; ++j) {
              u_rjk(r,j,k) += u_rsk_value * legendreBaseValues[1][s][j];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t>& u_ijk = legendre;
    {
      for (size_t r=0; r < gaussLobatto[0]->numberOfPoints(); ++r) {
        const real_t& wr = gaussLobatto[0]->weight(r);
        for (size_t j=0; j < __legendreDimensions[1]; ++j) {
          for (size_t k=0; k < __legendreDimensions[2]; ++k) {
            const real_t u_rjk_value = u_rjk(r,j,k) * wr;
            for (size_t i=0; i<__legendreDimensions[0]; ++i) {
              u_ijk(i,j,k) += u_rjk_value * legendreBaseValues[0][r][i];
            }
          }
        }
      }
    }
    
    //    return;
    for (size_t i=0; i < __legendreDimensions[0]; ++i) {
      const real_t wi_8 = 1./8 * (2*i+1);
      for (size_t j=0;j < __legendreDimensions[1]; ++j) { 
        const real_t wi_wj_8 = wi_8 * (2*j+1);
        for (size_t k=0; k < __legendreDimensions[2]; ++k) {
          const real_t wi_wj_wk_8 = wi_wj_8 * (2*k+1);
          u_ijk(i,j,k) *= wi_wj_wk_8 ;
        }
      }
    }
  }

  void __lagrangeToLegendreTransposed(const Structured3DVector<real_t>& legendre,
                                      Structured3DVector<real_t>& lagrange) const
  {
    lagrange = 0;
    
    TinyVector<3, ConstReferenceCounting<LegendreBasis> > legendreBasis;
    TinyVector<3, ConstReferenceCounting<LagrangeBasis> > lagrangeBasis;
    for (size_t i=0; i<3; ++i) {
      legendreBasis[i] = new LegendreBasis(__legendreDegrees[i]);
      lagrangeBasis[i] = new LagrangeBasis(__lagrangeVertices[i]);
    }
    
    TinyVector<3, ConstReferenceCounting<GaussLobatto> > gaussLobatto;
    for (size_t i=0; i<3; ++i) {
      gaussLobatto[i] = GaussLobattoManager::instance().getReference((__legendreDegrees[i]+__lagrangeDegrees[i]) + 1);
    }

    TinyVector<3,Vector<Vector<real_t> > >  legendreBaseValues;
    TinyVector<3,Vector<Vector<real_t> > >  lagrangeBaseValues;
    for (size_t i=0; i<3; ++i) {
      legendreBaseValues[i].resize(gaussLobatto[i]->numberOfPoints());
      lagrangeBaseValues[i].resize(gaussLobatto[i]->numberOfPoints());
    }
    
    for(size_t i=0; i<3; ++i) {
      const GaussLobatto& gaussLobattoPoints = *gaussLobatto[i];
      for(size_t j=0; j<gaussLobattoPoints.numberOfPoints(); ++j){
        
        // Compute position of vertex j in direction i
        const real_t& xj = gaussLobattoPoints(j);
        
        legendreBaseValues[i][j].resize(legendreBasis[i]->dimension());
        legendreBasis[i]->getValues(xj,legendreBaseValues[i][j]);
        
        lagrangeBaseValues[i][j].resize(lagrangeBasis[i]->dimension());
        lagrangeBasis[i]->getValues(xj,lagrangeBaseValues[i][j]);
      }
    }

    Structured3DVector<real_t> u_ijk = legendre;

    for (size_t i=0; i < __legendreDimensions[0]; ++i) {
      const real_t wi_8 = 1./8 * (2*i+1);
      for (size_t j=0;j < __legendreDimensions[1]; ++j) { 
        const real_t wi_wj_8 = wi_8 * (2*j+1);
        for (size_t k=0; k < __legendreDimensions[2]; ++k) {
          const real_t wi_wj_wk_8 = wi_wj_8 * (2*k+1);
          u_ijk(i,j,k) *= wi_wj_wk_8 ;
        }
      }
    }

    
    Structured3DVector<real_t> u_ijt(Array3DShape(__legendreDimensions[0],
                                                  __legendreDimensions[1],
                                                  gaussLobatto[2]->numberOfPoints()));
    u_ijt = 0;

    {
      for (size_t i=0; i < __legendreDimensions[0]; ++i) {
        for (size_t j=0; j < __legendreDimensions[1]; ++j) {
          for (size_t k=0; k < __legendreDimensions[2]; ++k) {
            const real_t& u_ijk_value = u_ijk(i,j,k);
            for (size_t t=0; t< gaussLobatto[2]->numberOfPoints(); ++t) {
              u_ijt(i,j,t) += u_ijk_value * legendreBaseValues[2][t][k];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_ist(Array3DShape(__legendreDimensions[0],
                                                   gaussLobatto[1]->numberOfPoints(),
                                                   gaussLobatto[2]->numberOfPoints()));
    u_ist = 0;
    
    {
      for (size_t i=0; i <__legendreDimensions[0]; ++i) {
        for (size_t j=0; j <__legendreDimensions[1]; ++j) {
          for (size_t t=0; t < gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& u_ijt_value = u_ijt(i,j,t) ;
            for (size_t s=0; s < gaussLobatto[1]->numberOfPoints(); ++s) {
              u_ist(i,s,t) += u_ijt_value * legendreBaseValues[1][s][j];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_rst(Array3DShape(gaussLobatto[0]->numberOfPoints(),
                                                  gaussLobatto[1]->numberOfPoints(),
                                                  gaussLobatto[2]->numberOfPoints()));
    u_rst = 0;
    {
      for (size_t i=0; i < __legendreDimensions[0]; ++i) {
        for (size_t s=0; s <gaussLobatto[1]->numberOfPoints(); ++s) {
          for (size_t t=0; t <gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& u_ist_value = u_ist(i,s,t);
            for (size_t r=0; r<gaussLobatto[0]->numberOfPoints(); ++r) {
              u_rst(r,s,t) += u_ist_value * legendreBaseValues[0][r][i];
            }
          }
        }
      }
    }

    Structured3DVector<real_t> u_rst0(Array3DShape(gaussLobatto[0]->numberOfPoints(),
                                                   gaussLobatto[1]->numberOfPoints(),
                                                   __lagrangeDimensions[2]));
    u_rst0 = 0;
    
    {
      for (size_t r=0; r < gaussLobatto[0]->numberOfPoints(); ++r) {
        for (size_t s=0; s < gaussLobatto[1]->numberOfPoints(); ++s) {
          for (size_t t=0; t < gaussLobatto[2]->numberOfPoints(); ++t) {
            const real_t& wt = gaussLobatto[2]->weight(t);
            const real_t u_rst_value = u_rst(r,s,t) * wt;
            for (size_t t0=0; t0<__lagrangeDimensions[2]; ++t0) {
              u_rst0(r,s,t0) += u_rst_value * lagrangeBaseValues[2][t][t0];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t> u_rs0t0(Array3DShape(gaussLobatto[0]->numberOfPoints(),
                                                    __lagrangeDimensions[1],
                                                    __lagrangeDimensions[2]));
    u_rs0t0 = 0;
    
    {
      for (size_t r=0; r < gaussLobatto[0]->numberOfPoints(); ++r) {
        for (size_t s=0; s < gaussLobatto[1]->numberOfPoints(); ++s) {
          const real_t& ws = gaussLobatto[1]->weight(s);
          for (size_t t0=0; t0 < __lagrangeDimensions[2]; ++t0) {
            const real_t u_rst0_value = u_rst0(r,s,t0) * ws;
            for (size_t s0=0; s0<__lagrangeDimensions[1]; ++s0) {
              u_rs0t0(r,s0,t0) += u_rst0_value * lagrangeBaseValues[1][s][s0];
            }
          }
        }
      }
    }
    
    Structured3DVector<real_t>& u_r0s0t0 = lagrange;
    {
      for (size_t r=0; r < gaussLobatto[0]->numberOfPoints(); ++r) {
        const real_t& wr = gaussLobatto[0]->weight(r);
        for (size_t s0=0; s0 < __lagrangeDimensions[1]; ++s0) {
          for (size_t t0=0; t0 < __lagrangeDimensions[2]; ++t0) {
            const real_t u_rs0t0_value = u_rs0t0(r,s0,t0) * wr;
            for (size_t r0=0; r0<__lagrangeDimensions[0]; ++r0) {
              u_r0s0t0(r0,s0,t0) += u_rs0t0_value * lagrangeBaseValues[0][r][r0];
            }
          }
        }
      }
    }
  }
  
  void __residuToLegendre (const Structured3DVector<real_t>& residu,
                           Structured3DVector<real_t>& legendre) const 
  {
    legendre = residu;

    for (size_t i=0; i < __legendreDimensions[0]; ++i) {
      const real_t wi_8 = 1./8 * (2*i+1);
      for (size_t j=0;j < __legendreDimensions[1]; ++j) { 
        const real_t wi_wj_8 = wi_8 * (2*j+1);
        for (size_t k=0; k < __legendreDimensions[2]; ++k) {
          const real_t wi_wj_wk_8 = wi_wj_8 * (2*k+1);
          legendre(i,j,k) *= wi_wj_wk_8 ;
        }
      }
    }
  }

  void __legendreToResidu(const Structured3DVector<real_t>& legendre,
                          Structured3DVector<real_t>& residu) const 
  {
    residu = legendre;

    for (size_t i=0; i < __legendreDimensions[0]; ++i) {
      const real_t wi_8 = 1./8 * (2*i+1);
      for (size_t j=0;j < __legendreDimensions[1]; ++j) { 
        const real_t wi_wj_8 = wi_8 * (2*j+1);
        for (size_t k=0; k < __legendreDimensions[2]; ++k) {
          const real_t wi_wj_wk_8 = wi_wj_8 * (2*k+1);
          residu(i,j,k) /= wi_wj_wk_8 ;
        }
      }
    }
  }

public:
  Internal(const Problem& problem)
    : __problem(problem)
  {
    ;
  }
  
  void computes(const Vector<real_t>& r,
                Vector<real_t>& z) const
  {
    const Mesh& solverMesh = SolverInformationCenter::instance().mesh();
    const DiscretizationType& discretization = SolverInformationCenter::instance().discretizationType();
    
    // Only works for scalar
    const ScalarDiscretizationTypeBase& scalarDiscretization = discretization[0];
    if (scalarDiscretization.type() != ScalarDiscretizationTypeBase::spectralLegendre) {
      throw ErrorHandler(__FILE__,__LINE__,
                         "bad discretization",
                         ErrorHandler::unexpected);
    }
    
    const ScalarDiscretizationTypeSpectral& spectralDiscretization
      = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(scalarDiscretization);
    
    __legendreDegrees = spectralDiscretization.degrees();
    __legendreDimensions = __legendreDegrees + TinyVector<3,size_t>(1,1,1);
    __a = spectralDiscretization.a();
    __b = spectralDiscretization.b();

    const bool useGaussLobatto = false;
    if(useGaussLobatto) {
      __lagrangeDegrees    = __legendreDegrees;
      __lagrangeDimensions = __legendreDimensions;

      for (size_t i=0; i<3; ++i) {
        __lagrangeVertices[i].resize(__legendreDimensions[i]);
      }

      for (size_t i=0; i<3; ++i) {
        const GaussLobatto& gaussLobatto = GaussLobattoManager::instance().get(__legendreDegrees[i]);
        __lagrangeVertices[i] = gaussLobatto.vertices();
      }
    } else {
      __lagrangeDegrees    = __legendreDegrees + TinyVector<3,size_t>(0,0,0);
      __lagrangeDimensions = __lagrangeDegrees + TinyVector<3,size_t>(1,1,1);

      for (size_t i=0; i<3; ++i) {
        __lagrangeVertices[i].resize(__lagrangeDimensions[i]);
      }

      for (size_t i=0; i<3; ++i) {

        const real_t h = std::abs(__b[i]-__a[i])/(__lagrangeDegrees[i]);
        Vector<real_t>& vertices = __lagrangeVertices[i];
        const real_t x0 = std::min(__a[i],__b[i]);

        for (size_t j=0; j<vertices.size(); ++j) {
          vertices[j] = x0 + j*h;
        }
      }
    }


    if (solverMesh.type() == Mesh::spectralMesh) {
      Structured3DVector<real_t> residu(__legendreDimensions);
      residu = r;
//       std::cout << "----------------------\n";
//       std::cout << "residu=\n" << residu;

      Structured3DVector<real_t> legendreResidu(__legendreDimensions);

      __residuToLegendre(residu,legendreResidu);
      //      legendreResidu = residu;
//       std::cout << "----------------------\n";
//       std::cout << "legendre=\n" << legendreResidu;

      Structured3DVector<real_t> lagrangeResidu(__lagrangeDimensions);

      __legendreToLagrange(legendreResidu, lagrangeResidu);

//       std::cout << "----------------------\n";
//       std::cout << "lagrange=\n" << lagrangeResidu;

      if (useGaussLobatto) {
        Structured3DMeshShape meshShape(__lagrangeDimensions-TinyVector<3,size_t>(2,2,2),
                                        __a,__b);
        ReferenceCounting<SpectralMesh> mesh
          = new SpectralMesh(meshShape,
                             new VerticesCorrespondance(__lagrangeDimensions[0]*__lagrangeDimensions[1]*__lagrangeDimensions[2]));
        __solveFEM<SpectralMesh>(mesh,lagrangeResidu);
      } else {
        Structured3DMeshShape meshShape(__lagrangeDimensions,
                                        __a,__b);
        ReferenceCounting<Structured3DMesh> mesh
          = new Structured3DMesh(meshShape,
                                 new VerticesCorrespondance(meshShape.numberOfVertices()));
        __solveFEM<Structured3DMesh>(mesh,lagrangeResidu);
      }

//       std::cout << "----------------------\n";
//       std::cout << "lagrange=\n" << lagrangeResidu << '\n';

      __lagrangeToLegendre(lagrangeResidu, legendreResidu);
//       std::cout << "----------------------\n";
//       std::cout << "legendre=\n" << legendreResidu << '\n';
      
      __legendreToResidu(legendreResidu,residu);

      z = residu;
//       std::cout << "----------------------\n";
//       std::cout << "residu=\n" << residu << '\n';
//       //      std::exit(0);
//       std::cout << "(z,r) 1 = " << r*z << '\n';
//       std::cout << Norm(z-r)/Norm(r) << '\n';
    } else {
      fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
      std::exit(1);
    }

  }

  void computesTransposed(const Vector<real_t>& r,
                          Vector<real_t>& z) const
  {
    const Mesh& solverMesh = SolverInformationCenter::instance().mesh();
    const DiscretizationType& discretization = SolverInformationCenter::instance().discretizationType();
    
    // Only works for scalar
    const ScalarDiscretizationTypeBase& scalarDiscretization = discretization[0];
    if (scalarDiscretization.type() != ScalarDiscretizationTypeBase::spectralLegendre) {
      throw ErrorHandler(__FILE__,__LINE__,
                         "bad discretization",
                         ErrorHandler::unexpected);
    }
    
    const ScalarDiscretizationTypeSpectral& spectralDiscretization
      = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(scalarDiscretization);
    
    __legendreDegrees = spectralDiscretization.degrees();
    __legendreDimensions = __legendreDegrees + TinyVector<3,size_t>(1,1,1);
    __a = spectralDiscretization.a();
    __b = spectralDiscretization.b();

    const bool useGaussLobatto = false;
    if(useGaussLobatto) {
      __lagrangeDegrees    = __legendreDegrees;
      __lagrangeDimensions = __legendreDimensions;

      for (size_t i=0; i<3; ++i) {
        __lagrangeVertices[i].resize(__legendreDimensions[i]);
      }

      for (size_t i=0; i<3; ++i) {
        const GaussLobatto& gaussLobatto = GaussLobattoManager::instance().get(__legendreDegrees[i]);
        __lagrangeVertices[i] = gaussLobatto.vertices();
      }
    } else {
      __lagrangeDegrees    = __legendreDegrees + TinyVector<3,size_t>(0,0,0);
      __lagrangeDimensions = __lagrangeDegrees + TinyVector<3,size_t>(1,1,1);

      for (size_t i=0; i<3; ++i) {
        __lagrangeVertices[i].resize(__lagrangeDimensions[i]);
      }

      for (size_t i=0; i<3; ++i) {

        const real_t h = std::abs(__b[i]-__a[i])/(__lagrangeDegrees[i]);
        Vector<real_t>& vertices = __lagrangeVertices[i];
        const real_t x0 = std::min(__a[i],__b[i]);

        for (size_t j=0; j<vertices.size(); ++j) {
          vertices[j] = x0 + j*h;
        }
      }
    }

    if (solverMesh.type() == Mesh::spectralMesh) {
      Structured3DVector<real_t> residu(__legendreDimensions);
      residu = r;

      Structured3DVector<real_t> legendreResidu(__legendreDimensions);
      Structured3DVector<real_t> lagrangeResidu(__lagrangeDimensions);

      __legendreToResidu(residu,legendreResidu);
      __lagrangeToLegendreTransposed(legendreResidu,lagrangeResidu);
      Structured3DMeshShape meshShape(__lagrangeDimensions,
                                      __a,__b);
      ReferenceCounting<Structured3DMesh> mesh
        = new Structured3DMesh(meshShape,
                               new VerticesCorrespondance(meshShape.numberOfVertices()));
      __solveFEM<Structured3DMesh>(mesh,lagrangeResidu);
      __legendreToLagrangeTransposed(lagrangeResidu, legendreResidu);
      __residuToLegendre(legendreResidu,residu);

      z = residu;
//       std::cout << "(z,r) 2 = " << r*z << '\n';
//       std::cout << Norm(z-r)/Norm(r) << '\n';
    } else {
      fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
      std::exit(1);
    }

  }
};

void SpectralFEMPreconditioner::
initializes()
{
  int oldVerbosity = StreamCenter::instance().getDebugLevel();
  StreamCenter::instance().setDebugLevel(oldVerbosity-5);
  __internal = new SpectralFEMPreconditioner::Internal(__problem);
  StreamCenter::instance().setDebugLevel(oldVerbosity);
}

void
SpectralFEMPreconditioner::
computes(const Vector<real_t>& r,
         Vector<real_t>& z) const
{
  int oldVerbosity = StreamCenter::instance().getDebugLevel();
  StreamCenter::instance().setDebugLevel(oldVerbosity-5);
  __internal->computes(r,z);
  StreamCenter::instance().setDebugLevel(oldVerbosity);
}

void
SpectralFEMPreconditioner::
computesTransposed(const Vector<real_t>& r,
         Vector<real_t>& z) const
{
  int oldVerbosity = StreamCenter::instance().getDebugLevel();
  StreamCenter::instance().setDebugLevel(oldVerbosity-5);
  __internal->computesTransposed(r,z);
  StreamCenter::instance().setDebugLevel(oldVerbosity);
}

SpectralFEMPreconditioner::
SpectralFEMPreconditioner(const Problem& problem)
  : Preconditioner (problem,
                    Preconditioner::spectralFEM)
{
  ;
}

SpectralFEMPreconditioner::
~SpectralFEMPreconditioner()
{
  ;
}

---