/*********************************************************************/
/* File:   vefc.cpp                                                  */
/* Author: Joachim Schoeberl                                         */
/* Date:   14. July. 2002                                            */
/*********************************************************************/

/* 
   VEFC high order Preconditioner
*/

#include <multigrid.hpp>

namespace ngmg
{
  using namespace ngmg;
  VEFC_Matrix :: 
  VEFC_Matrix (const BaseMatrix & amat, 
	       const FESpace & afes)
    : mat(amat), fes(afes)
  {
    cout << "update VEFC Matrix, fes = " << typeid(fes).name() << endl;
    Table<int> * blocks = fes.CreateSmoothingBlocks();

    jacobi = dynamic_cast<const BaseSparseMatrix&> (mat)
      .CreateBlockJacobiPrecond(*blocks);
  }
  
  VEFC_Matrix :: ~VEFC_Matrix ()
  {
    ;
  }


  void VEFC_Matrix :: Mult (const BaseVector & x, BaseVector & y) const
  {
    int i;

    BaseVector & hy = *mat.CreateVector();

    ARRAY<int> dofs;
    int ne = fes.GetMeshAccess().GetNE();
    
    LocalHeap lh(1000);

    y = x;
    hy = 0;
    for (i = 0; i < ne; i++)
      {
	fes.GetDofNrs (i, dofs);

	Vector<> elv(dofs.Size());
	y.GetIndirect (dofs, elv);

	const FE_Augmented_TetP * tet = 
	  dynamic_cast<const FE_Augmented_TetP*> (&fes.GetFE(i, lh));
	if (tet) tet->UnifyTrans (elv);

	const FE_Augmented_TrigP * trig = 
	  dynamic_cast<const FE_Augmented_TrigP*> (&fes.GetFE(i, lh));
	if (trig) trig->UnifyTrans (elv);
	
	lh.CleanUp();

	/*
	try {
	  dynamic_cast<const FE_Augmented_TetP&> (fes.GetFE(i)).UnifyTrans (elv);
	}
	catch (...) { ; }

	try {
	  dynamic_cast<const FE_Augmented_TrigP&> (fes.GetFE(i)).UnifyTrans (elv);
	}
	catch (...) { ; }
	*/

	
	hy.AddIndirect (dofs, elv);
	elv = 0;
	y.SetIndirect (dofs, elv);
      }

    y = 0;
    jacobi->GSSmooth (y, hy);
    jacobi->GSSmoothBack (y, hy);

    //    jacobi->Mult (x, y);

    
    hy = y;
    for (i = 0; i < ne; i++)
      {
	fes.GetDofNrs (i, dofs);
	Vector<> elv(dofs.Size());
	hy.GetIndirect (dofs, elv);

	const FE_Augmented_TetP * tet = 
	  dynamic_cast<const FE_Augmented_TetP*> (&fes.GetFE(i, lh));
	if (tet) tet->Unify (elv);

	const FE_Augmented_TrigP * trig = 
	  dynamic_cast<const FE_Augmented_TrigP*> (&fes.GetFE(i, lh));
	if (trig) trig->Unify (elv);

	lh.CleanUp();
	/*	
	try {
	  dynamic_cast<const FE_Augmented_TetP&> (fes.GetFE(i)).Unify (elv);
	}
	catch (...) { ; }

	try {
	  dynamic_cast<const FE_Augmented_TrigP&> (fes.GetFE(i)).Unify (elv);
	}
	catch (...) { ; }
	*/
	y.SetIndirect (dofs, elv);
      }
    

    delete &hy;
    
    
    //    jacobi->Mult (x, y);

    /*
    BaseVector & cres = *cpre->CreateVector();
    BaseVector & cw = *cpre->CreateVector();
    BaseVector & res = *CreateVector();

    y = 0;
    //  jacsmoother->GSSmooth (y, x);
    smoother->PreSmooth (level, y, x, smoothingsteps);

    res = x - (*mat) * y;  // mat->Residuum (y, x, res); 
    cres = res.Range (0, cres.Size()); // res.GetPart (1, cres);
    cw = (*cpre) * cres;
    y.Range (0, cw.Size()) += cw;

    smoother->PostSmooth (level, y, x, smoothingsteps);
    // jacsmoother->GSSmoothBack (y, x);

    delete &cres;
    delete &cw;
    delete &res;
    */
  }

  ngla::BaseVector * VEFC_Matrix :: CreateVector () const
  {
    return mat.CreateVector();
  }

  ostream & VEFC_Matrix :: Print (ostream & s) const
  {
    s << "Twolevel Preconditioner\n";
    return s;
  }




  void VEFC_Matrix :: MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    //    if (cpre) cpre->MemoryUsage (mu);
  }


}