SpectralFEMPreconditioner::Internal Class Reference

Collaboration diagram for SpectralFEMPreconditioner::Internal:

[legend]

List of all members.

Public Member Functions
	Internal (const Problem &problem)
void	computes (const Vector< real_t > &r, Vector< real_t > &z) const
void	computesTransposed (const Vector< real_t > &r, Vector< real_t > &z) const
Private Member Functions
template<typename MeshType>
void	__solveFEM (ConstReferenceCounting< MeshType > mesh, Structured3DVector< real_t > &lagrange) const
void	__legendreToLagrange (const Structured3DVector< real_t > &legendre, Structured3DVector< real_t > &lagrange) const
void	__legendreToLagrangeTransposed (const Structured3DVector< real_t > &lagrange, Structured3DVector< real_t > &legendre) const
void	__lagrangeToLegendre (const Structured3DVector< real_t > &lagrange, Structured3DVector< real_t > &legendre) const
void	__lagrangeToLegendreTransposed (const Structured3DVector< real_t > &legendre, Structured3DVector< real_t > &lagrange) const
void	__residuToLegendre (const Structured3DVector< real_t > &residu, Structured3DVector< real_t > &legendre) const
void	__legendreToResidu (const Structured3DVector< real_t > &legendre, Structured3DVector< real_t > &residu) const
Private Attributes
const Problem &	__problem
TinyVector< 3, size_t >	__legendreDegrees
TinyVector< 3, size_t >	__legendreDimensions
TinyVector< 3, size_t >	__lagrangeDegrees
TinyVector< 3, size_t >	__lagrangeDimensions
TinyVector< 3, real_t >	__a
TinyVector< 3, real_t >	__b
TinyVector< 3, Vector< real_t > >	__lagrangeVertices

---

Detailed Description

Definition at line 60 of file SpectralFEMPreconditioner.cpp.

---

Constructor & Destructor Documentation

SpectralFEMPreconditioner::Internal::Internal

(

const Problem &

problem

)

[inline]

Constructor

Parameters:

problem

problem to precondition

Definition at line 708 of file SpectralFEMPreconditioner.cpp.

    : __problem(problem)
  {
    ;
  }

---

Member Function Documentation

template<typename MeshType>

void SpectralFEMPreconditioner::Internal::__solveFEM	(	ConstReferenceCounting< MeshType >	mesh,
		Structured3DVector< real_t > &	lagrange
	)			const [inline, private]

Definition at line 77 of file SpectralFEMPreconditioner.cpp.

References __problem, DiscretizationType::add(), VariationalProblem::add(), FEMFunctionBuilder::build(), MultiGrid::computes(), DegreeOfFreedomSetBuilder::degreeOfFreedomSet(), BaseMatrix::doubleHashedMatrix, ffout(), FEMFunctionBuilder::getBuiltScalarFunction(), StaticBase< ParameterCenter >::instance(), ScalarDiscretizationTypeBase::lagrangianFEM1, VariationalOperator::normal, VariationalProblem::numberOfUnknown(), MemoryManager::ReserveMatrix(), ParameterCenter::set(), DegreeOfFreedomSet::size(), Timer::start(), Timer::stop(), Problem::type(), ErrorHandler::unexpected, and Problem::variational.

  {
    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"; 
  }

Here is the call graph for this function:

void SpectralFEMPreconditioner::Internal::__legendreToLagrange	(	const Structured3DVector< real_t > &	legendre,
		Structured3DVector< real_t > &	lagrange
	)			const [inline, private]

Definition at line 190 of file SpectralFEMPreconditioner.cpp.

References __lagrangeDimensions, __lagrangeVertices, __legendreDegrees, __legendreDimensions, and Vector< T >::size().

Referenced by computes().

  {
    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];
          }
        }
      }
    }
  }

Here is the call graph for this function:

void SpectralFEMPreconditioner::Internal::__legendreToLagrangeTransposed	(	const Structured3DVector< real_t > &	lagrange,
		Structured3DVector< real_t > &	legendre
	)			const [inline, private]

Definition at line 267 of file SpectralFEMPreconditioner.cpp.

References __lagrangeDimensions, __lagrangeVertices, __legendreDegrees, __legendreDimensions, and Vector< T >::size().

Referenced by computesTransposed().

  {
    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];
          }
        }
      }
    }
  }

Here is the call graph for this function:

void SpectralFEMPreconditioner::Internal::__lagrangeToLegendre	(	const Structured3DVector< real_t > &	lagrange,
		Structured3DVector< real_t > &	legendre
	)			const [inline, private]

Taking into account of the renormalization of the legendre polynomials

Definition at line 344 of file SpectralFEMPreconditioner.cpp.

References __lagrangeDegrees, __lagrangeDimensions, __lagrangeVertices, __legendreDegrees, __legendreDimensions, GaussLobattoManager::getReference(), StaticBase< GaussLobattoManager >::instance(), and GaussLobatto::numberOfPoints().

Referenced by computes().

  {
    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 ;
        }
      }
    }
  }

Here is the call graph for this function:

void SpectralFEMPreconditioner::Internal::__lagrangeToLegendreTransposed	(	const Structured3DVector< real_t > &	legendre,
		Structured3DVector< real_t > &	lagrange
	)			const [inline, private]

Taking into account of the renormalization of the legendre polynomials

Definition at line 505 of file SpectralFEMPreconditioner.cpp.

References __lagrangeDegrees, __lagrangeDimensions, __lagrangeVertices, __legendreDegrees, __legendreDimensions, GaussLobattoManager::getReference(), StaticBase< GaussLobattoManager >::instance(), and GaussLobatto::numberOfPoints().

Referenced by computesTransposed().

  {
    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];
            }
          }
        }
      }
    }
  }

Here is the call graph for this function:

void SpectralFEMPreconditioner::Internal::__residuToLegendre	(	const Structured3DVector< real_t > &	residu,
		Structured3DVector< real_t > &	legendre
	)			const [inline, private]

Definition at line 668 of file SpectralFEMPreconditioner.cpp.

References __legendreDimensions.

Referenced by computes(), and computesTransposed().

  {
    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 SpectralFEMPreconditioner::Internal::__legendreToResidu	(	const Structured3DVector< real_t > &	legendre,
		Structured3DVector< real_t > &	residu
	)			const [inline, private]

Definition at line 685 of file SpectralFEMPreconditioner.cpp.

References __legendreDimensions.

Referenced by computes(), and computesTransposed().

  {
    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 ;
        }
      }
    }
  }

void SpectralFEMPreconditioner::Internal::computes	(	const Vector< real_t > &	r,
		Vector< real_t > &	z
	)			const [inline]

Computes

Parameters:

	r	given residue
	z	preconditioned residue

Definition at line 720 of file SpectralFEMPreconditioner.cpp.

References __a, __b, __lagrangeDegrees, __lagrangeDimensions, __lagrangeToLegendre(), __lagrangeVertices, __legendreDegrees, __legendreDimensions, __legendreToLagrange(), __legendreToResidu(), __residuToLegendre(), ScalarDiscretizationTypeSpectral::a(), ScalarDiscretizationTypeSpectral::b(), ScalarDiscretizationTypeSpectral::degrees(), SolverInformationCenter::discretizationType(), fferr(), GaussLobattoManager::get(), StaticBase< GaussLobattoManager >::instance(), StaticBase< SolverInformationCenter >::instance(), SolverInformationCenter::mesh(), Structured3DMeshShape::numberOfVertices(), Vector< T >::size(), ScalarDiscretizationTypeBase::spectralLegendre, Mesh::spectralMesh, Mesh::type(), ScalarDiscretizationTypeBase::type(), ErrorHandler::unexpected, and GaussLobatto::vertices().

  {
    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);
    }

  }

Here is the call graph for this function:

void SpectralFEMPreconditioner::Internal::computesTransposed	(	const Vector< real_t > &	r,
		Vector< real_t > &	z
	)			const [inline]

Definition at line 834 of file SpectralFEMPreconditioner.cpp.

References __a, __b, __lagrangeDegrees, __lagrangeDimensions, __lagrangeToLegendreTransposed(), __lagrangeVertices, __legendreDegrees, __legendreDimensions, __legendreToLagrangeTransposed(), __legendreToResidu(), __residuToLegendre(), ScalarDiscretizationTypeSpectral::a(), ScalarDiscretizationTypeSpectral::b(), ScalarDiscretizationTypeSpectral::degrees(), SolverInformationCenter::discretizationType(), fferr(), GaussLobattoManager::get(), StaticBase< GaussLobattoManager >::instance(), StaticBase< SolverInformationCenter >::instance(), SolverInformationCenter::mesh(), Structured3DMeshShape::numberOfVertices(), Vector< T >::size(), ScalarDiscretizationTypeBase::spectralLegendre, Mesh::spectralMesh, Mesh::type(), ScalarDiscretizationTypeBase::type(), ErrorHandler::unexpected, and GaussLobatto::vertices().

  {
    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);
    }

  }

Here is the call graph for this function:

---

Member Data Documentation

const Problem& SpectralFEMPreconditioner::Internal::__problem [private]

Definition at line 63 of file SpectralFEMPreconditioner.cpp.

Referenced by __solveFEM().

TinyVector<3,size_t> SpectralFEMPreconditioner::Internal::__legendreDegrees [mutable, private]

Definition at line 65 of file SpectralFEMPreconditioner.cpp.

Referenced by __lagrangeToLegendre(), __lagrangeToLegendreTransposed(), __legendreToLagrange(), __legendreToLagrangeTransposed(), computes(), and computesTransposed().

TinyVector<3,size_t> SpectralFEMPreconditioner::Internal::__legendreDimensions [mutable, private]

Definition at line 66 of file SpectralFEMPreconditioner.cpp.

Referenced by __lagrangeToLegendre(), __lagrangeToLegendreTransposed(), __legendreToLagrange(), __legendreToLagrangeTransposed(), __legendreToResidu(), __residuToLegendre(), computes(), and computesTransposed().

TinyVector<3,size_t> SpectralFEMPreconditioner::Internal::__lagrangeDegrees [mutable, private]

Definition at line 68 of file SpectralFEMPreconditioner.cpp.

Referenced by __lagrangeToLegendre(), __lagrangeToLegendreTransposed(), computes(), and computesTransposed().

TinyVector<3,size_t> SpectralFEMPreconditioner::Internal::__lagrangeDimensions [mutable, private]

Definition at line 69 of file SpectralFEMPreconditioner.cpp.

Referenced by __lagrangeToLegendre(), __lagrangeToLegendreTransposed(), __legendreToLagrange(), __legendreToLagrangeTransposed(), computes(), and computesTransposed().

TinyVector<3,real_t> SpectralFEMPreconditioner::Internal::__a [mutable, private]

Definition at line 71 of file SpectralFEMPreconditioner.cpp.

Referenced by computes(), and computesTransposed().

TinyVector<3,real_t> SpectralFEMPreconditioner::Internal::__b [mutable, private]

Definition at line 72 of file SpectralFEMPreconditioner.cpp.

Referenced by computes(), and computesTransposed().

TinyVector<3,Vector<real_t> > SpectralFEMPreconditioner::Internal::__lagrangeVertices [mutable, private]

Definition at line 74 of file SpectralFEMPreconditioner.cpp.

Referenced by __lagrangeToLegendre(), __lagrangeToLegendreTransposed(), __legendreToLagrange(), __legendreToLagrangeTransposed(), computes(), and computesTransposed().

---

The documentation for this class was generated from the following file:

• SpectralFEMPreconditioner.cpp