MultiGrid Class Reference

#include <MultiGrid.hpp>

Inheritance diagram for MultiGrid:

[legend]

Collaboration diagram for MultiGrid:

[legend]

List of all members.

Public Types
enum	Type {   none, diagonale, incompleteCholeskiFactorization, multigrid,   spectralFEM }
Public Member Functions
std::string	name () const
	MultiGrid (const Problem &problem, const SparseMatrix &A, const DegreeOfFreedomSet &degreeOfFreedomSet)
void	copyDirichlet (const Vector< real_t > &u, Vector< real_t > &v, const size_t level) const
void	residuDirichlet (Vector< real_t > &v, const Vector< real_t > &u1, const Vector< real_t > &u2, const size_t level) const
void	setDirichlet (const size_t level, Vector< real_t > &v) const
void	getDirichlet (const size_t level, Vector< real_t > &v) const
void	interpolateDirichlet (Vector< real_t > &u2, const Vector< real_t > &u1, const size_t currentLevel) const
void	restrictDirichlet (const Vector< real_t > &u2, Vector< real_t > &u1, const size_t currentLevel) const
void	initializes ()
	Initialization of the Preconditioner.
void	weightJacobi (Structured3DVector< real_t > &u, const Structured3DVector< real_t > &b, const size_t currentLevel) const
void	computeResidu (Structured3DVector< real_t > &r, const Structured3DVector< real_t > &u, const Structured3DVector< real_t > &f, const size_t currentLevel) const
void	correct (Structured3DVector< real_t > &u, const Structured3DVector< real_t > &u2, const size_t currentLevel) const
void	project (Structured3DVector< real_t > &u2, const Structured3DVector< real_t > &u, const size_t currentLevel) const
void	multigrid (Structured3DVector< real_t > &u, const Structured3DVector< real_t > &b, int currentLevel) const
void	FMV (Structured3DVector< real_t > &u, const Structured3DVector< real_t > &b, int currentLevel) const
void	computes (const Vector< real_t > &f, Vector< real_t > &v) const
	Computes .
const Preconditioner::Type &	Type () const
virtual void	computesTransposed (const Vector< real_t > &r, Vector< real_t > &z) const
	Computes .
Protected Attributes
Preconditioner::Type	__type
	The type of preconditioner.
const Problem &	__problem
Private Attributes
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
std::vector< Vector< real_t > >	__dirichletValues
Classes
class	DirichletPositions

---

Detailed Description

Multi Grid method.

Resolve the linear system #0.

Author:

Stéphane Del Pino.

Definition at line 60 of file MultiGrid.hpp.

---

Member Enumeration Documentation

enum Preconditioner::Type [inherited]

Enumerator:

none
diagonale
incompleteCholeskiFactorization
multigrid
spectralFEM

Definition at line 48 of file Preconditioner.hpp.

            {
    none,
    diagonale,
    incompleteCholeskiFactorization,
    multigrid,
    spectralFEM
  };

---

Constructor & Destructor Documentation

MultiGrid::MultiGrid	(	const Problem &	problem,
		const SparseMatrix &	A,
		const DegreeOfFreedomSet &	degreeOfFreedomSet
	)			[inline]

Definition at line 142 of file MultiGrid.hpp.

References __degreeOfFreedomSet, __dirichlet, __dirichletValues, __level, __meshesShape, __meshShape, __options, __positionList, Preconditioner::__problem, Structured3DMeshShape::a(), Structured3DMeshShape::b(), Mesh::cartesianHexahedraMesh, SolverInformationCenter::discretizationType(), MultiGridOptions::epsilon(), epsilon, Structured3DMeshShape::hx(), Structured3DMeshShape::hy(), Structured3DMeshShape::hz(), StaticBase< SolverInformationCenter >::instance(), ScalarDiscretizationTypeBase::lagrangianFEM1, MultiGridOptions::maxiter(), maxiter, SolverInformationCenter::mesh(), MultiGridOptions::mu1(), mu1, MultiGridOptions::mu2(), mu2, ErrorHandler::normal, MultiGridOptions::nu1(), nu1, MultiGridOptions::nu2(), nu2, Structured3DMeshShape::numberOfVertices(), Structured3DVector< T >::nx(), Structured3DMeshShape::nx(), Array3DShape::nx(), Structured3DVector< T >::ny(), Structured3DMeshShape::ny(), Array3DShape::ny(), Structured3DVector< T >::nz(), Structured3DMeshShape::nz(), Array3DShape::nz(), MultiGridOptions::omega(), omega, Structured3DVector< T >::resize(), MultiGrid::DirichletPositions::setPosition(), Structured3DVector< T >::shape(), SpectralMesh::shape(), Structured3DMeshShape::shape(), Structured3DMesh::shape(), DegreeOfFreedomSet::size(), Mesh::spectralMesh, Mesh::type(), and GetParameter< T >::value().

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

Here is the call graph for this function:

---

Member Function Documentation

std::string MultiGrid::name

(

)

const [inline, virtual]

Implements Preconditioner.

Definition at line 136 of file MultiGrid.hpp.

  {
    return "finite differences multigrid";
  }

void MultiGrid::copyDirichlet	(	const Vector< real_t > &	u,
		Vector< real_t > &	v,
		const size_t	level
	)			const [inline]

Definition at line 268 of file MultiGrid.hpp.

References __level, __meshesShape, __positionList, ASSERT, and Vector< T >::size().

Referenced by weightJacobi().

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

Here is the call graph for this function:

void MultiGrid::residuDirichlet	(	Vector< real_t > &	v,
		const Vector< real_t > &	u1,
		const Vector< real_t > &	u2,
		const size_t	level
	)			const [inline]

Definition at line 283 of file MultiGrid.hpp.

References __level, __meshesShape, __positionList, ASSERT, and Vector< T >::size().

Referenced by computeResidu().

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

Here is the call graph for this function:

void MultiGrid::setDirichlet	(	const size_t	level,
		Vector< real_t > &	v
	)			const [inline]

Definition at line 297 of file MultiGrid.hpp.

References __dirichletValues, __level, and __positionList.

Referenced by computes(), and correct().

  {
    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 MultiGrid::getDirichlet	(	const size_t	level,
		Vector< real_t > &	v
	)			const [inline]

Definition at line 307 of file MultiGrid.hpp.

References __dirichletValues, __level, and __positionList.

Referenced by computes(), and correct().

  {
    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 MultiGrid::interpolateDirichlet	(	Vector< real_t > &	u2,
		const Vector< real_t > &	u1,
		const size_t	currentLevel
	)			const [inline]

Definition at line 317 of file MultiGrid.hpp.

References __level, and __positionList.

Referenced by project().

  {
    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 MultiGrid::restrictDirichlet	(	const Vector< real_t > &	u2,
		Vector< real_t > &	u1,
		const size_t	currentLevel
	)			const [inline]

Definition at line 333 of file MultiGrid.hpp.

References __level, and __positionList.

  {
    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 MultiGrid::initializes

(

)

[inline, virtual]

Initialization of the Preconditioner.

Implements Preconditioner.

Definition at line 348 of file MultiGrid.hpp.

  {
    ;
  }

void MultiGrid::weightJacobi	(	Structured3DVector< real_t > &	u,
		const Structured3DVector< real_t > &	b,
		const size_t	currentLevel
	)			const [inline]

Definition at line 353 of file MultiGrid.hpp.

References __meshShape, copyDirichlet(), Structured3DMeshShape::dimension(), Array3DShape::nx(), Array3DShape::ny(), Array3DShape::nz(), omega, Structured3DVector< T >::shape(), and Vector< T >::size().

Referenced by multigrid().

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

Here is the call graph for this function:

void MultiGrid::computeResidu	(	Structured3DVector< real_t > &	r,
		const Structured3DVector< real_t > &	u,
		const Structured3DVector< real_t > &	f,
		const size_t	currentLevel
	)			const [inline]

RECOMPUTE THOSE VALUES !

Definition at line 388 of file MultiGrid.hpp.

References __meshShape, ASSERT, Structured3DMeshShape::dimension(), Structured3DVector< T >::nx(), Structured3DVector< T >::ny(), Structured3DVector< T >::nz(), residuDirichlet(), and Structured3DVector< T >::shape().

Referenced by computes(), FMV(), and multigrid().

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

Here is the call graph for this function:

void MultiGrid::correct	(	Structured3DVector< real_t > &	u,
		const Structured3DVector< real_t > &	u2,
		const size_t	currentLevel
	)			const [inline]

Definition at line 421 of file MultiGrid.hpp.

References ASSERT, getDirichlet(), Structured3DVector< T >::nx(), Structured3DVector< T >::ny(), Structured3DVector< T >::nz(), setDirichlet(), and sum().

Referenced by FMV(), and multigrid().

  {

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

Here is the call graph for this function:

void MultiGrid::project	(	Structured3DVector< real_t > &	u2,
		const Structured3DVector< real_t > &	u,
		const size_t	currentLevel
	)			const [inline]

Definition at line 464 of file MultiGrid.hpp.

References interpolateDirichlet(), Structured3DVector< T >::nx(), Structured3DVector< T >::ny(), and Structured3DVector< T >::nz().

Referenced by FMV(), and multigrid().

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

Here is the call graph for this function:

void MultiGrid::multigrid	(	Structured3DVector< real_t > &	u,
		const Structured3DVector< real_t > &	b,
		int	currentLevel
	)			const [inline]

proceed to some more relaxation for the coasest grid

Definition at line 500 of file MultiGrid.hpp.

References computeResidu(), correct(), mu1, nu1, nu2, Structured3DVector< T >::nx(), Structured3DVector< T >::ny(), Structured3DVector< T >::nz(), project(), Structured3DVector< T >::shape(), and weightJacobi().

Referenced by FMV().

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

Here is the call graph for this function:

void MultiGrid::FMV	(	Structured3DVector< real_t > &	u,
		const Structured3DVector< real_t > &	b,
		int	currentLevel
	)			const [inline]

Definition at line 537 of file MultiGrid.hpp.

References computeResidu(), correct(), mu2, multigrid(), Structured3DVector< T >::nx(), Structured3DVector< T >::ny(), Structured3DVector< T >::nz(), project(), and Structured3DVector< T >::shape().

Referenced by computes().

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

Here is the call graph for this function:

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

Computes .

Implements Preconditioner.

Definition at line 563 of file MultiGrid.hpp.

References __level, __meshShape, ASSERT, computeResidu(), epsilon, ffout(), FMV(), getDirichlet(), Structured3DMeshShape::hx(), Structured3DMeshShape::hy(), Structured3DMeshShape::hz(), maxiter, setDirichlet(), Structured3DVector< T >::shape(), Structured3DMeshShape::shape(), and Vector< T >::size().

Referenced by SpectralFEMPreconditioner::Internal::__solveFEM().

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

Here is the call graph for this function:

const Preconditioner::Type& Preconditioner::Type

(

)

const [inline, inherited]

Definition at line 100 of file Preconditioner.hpp.

References Preconditioner::__type.

  {
    return __type;
  }

virtual void Preconditioner::computesTransposed	(	const Vector< real_t > &	r,
		Vector< real_t > &	z
	)			const [inline, virtual, inherited]

Computes .

Reimplemented in SpectralFEMPreconditioner.

Definition at line 72 of file Preconditioner.hpp.

References Preconditioner::computes().

  {
    this->computes(r,z);
  }

Here is the call graph for this function:

---

Member Data Documentation

size_t MultiGrid::__level [mutable, private]

Definition at line 111 of file MultiGrid.hpp.

Referenced by computes(), copyDirichlet(), getDirichlet(), interpolateDirichlet(), MultiGrid(), residuDirichlet(), restrictDirichlet(), and setDirichlet().

int MultiGrid::nu1 [private]

Definition at line 112 of file MultiGrid.hpp.

Referenced by multigrid(), and MultiGrid().

int MultiGrid::nu2 [private]

Definition at line 113 of file MultiGrid.hpp.

Referenced by multigrid(), and MultiGrid().

int MultiGrid::mu1 [private]

Definition at line 114 of file MultiGrid.hpp.

Referenced by multigrid(), and MultiGrid().

int MultiGrid::mu2 [private]

Definition at line 115 of file MultiGrid.hpp.

Referenced by FMV(), and MultiGrid().

int MultiGrid::maxiter [private]

Definition at line 117 of file MultiGrid.hpp.

Referenced by computes(), and MultiGrid().

real_t MultiGrid::omega [private]

Definition at line 119 of file MultiGrid.hpp.

Referenced by MultiGrid(), and weightJacobi().

real_t MultiGrid::epsilon [private]

Definition at line 120 of file MultiGrid.hpp.

Referenced by computes(), and MultiGrid().

const DegreeOfFreedomSet& MultiGrid::__degreeOfFreedomSet [private]

Definition at line 122 of file MultiGrid.hpp.

Referenced by MultiGrid().

Structured3DMeshShape MultiGrid::__meshShape [private]

Definition at line 124 of file MultiGrid.hpp.

Referenced by computeResidu(), computes(), MultiGrid(), and weightJacobi().

GetParameter<MultiGridOptions> MultiGrid::__options [private]

Definition at line 126 of file MultiGrid.hpp.

Referenced by MultiGrid().

Structured3DVector<bool> MultiGrid::__dirichlet [private]

Definition at line 128 of file MultiGrid.hpp.

Referenced by MultiGrid().

std::vector<Structured3DMeshShape> MultiGrid::__meshesShape [private]

Definition at line 130 of file MultiGrid.hpp.

Referenced by copyDirichlet(), MultiGrid(), and residuDirichlet().

std::list<DirichletPositions> MultiGrid::__positionList [private]

Definition at line 131 of file MultiGrid.hpp.

Referenced by copyDirichlet(), getDirichlet(), interpolateDirichlet(), MultiGrid(), residuDirichlet(), restrictDirichlet(), and setDirichlet().

std::vector<Vector<real_t> > MultiGrid::__dirichletValues [mutable, private]

Definition at line 133 of file MultiGrid.hpp.

Referenced by getDirichlet(), MultiGrid(), and setDirichlet().

Preconditioner::Type Preconditioner::__type [protected, inherited]

The type of preconditioner.

Definition at line 58 of file Preconditioner.hpp.

Referenced by Preconditioner::Type().

const Problem& Preconditioner::__problem [protected, inherited]

Definition at line 61 of file Preconditioner.hpp.

Referenced by SpectralFEMPreconditioner::initializes(), and MultiGrid().

---

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

• MultiGrid.hpp