/*********************************************************************/
/* File:   hdivfes.cpp                                               */
/* Author: Joachim Schoeberl                                         */
/* Date:   27. Jan. 2003                                             */
/*********************************************************************/

/* 
   Finite Element Space for H(div)
*/

#include <comp.hpp>

namespace ngcomp
{
  using namespace ngcomp;

  RaviartThomasFESpace :: 
  RaviartThomasFESpace (const MeshAccess & ama,
			int adim, bool acomplex)
    
    : FESpace (ama, 1, adim, acomplex)
  {
    /*
    tet     = 
    prism   = 
    pyramid = 
    */
    trig    = new FE_RTTrig0;
    /*
    quad    = 
    segm    = 
    */
    
    // CreateVecObject1(prol, EdgeProlongation,
    // dimension, iscomplex, *this);


    if (ma.GetDimension() == 2)
      {
	ARRAY<CoefficientFunction*> coeffs(1);
	coeffs[1] = new ConstantCoefficientFunction(1);
	evaluator = GetIntegrators().CreateBFI("masshdiv", 2, coeffs);
      }
  }
  
  RaviartThomasFESpace :: ~RaviartThomasFESpace ()
  {
    ;
  }


  FESpace * RaviartThomasFESpace :: 
  Create (const MeshAccess & ma, const Flags & flags)
  {
    int order = int(flags.GetNumFlag ("order", 0));
    int dim = int(flags.GetNumFlag ("dim", 1));
    bool iscomplex = flags.GetDefineFlag ("complex");

    if (order == 0)
      return new RaviartThomasFESpace (ma, dim, iscomplex);      
    else
      return new HDivHighOrderFESpace (ma,flags);
  }



  
  void RaviartThomasFESpace :: Update()
  {
    const MeshAccess & ma = GetMeshAccess();
    int level = ma.GetNLevels();
    
    if (level == ndlevel.Size())
      return;
    
    if (ma.GetDimension() == 2)
      ndlevel.Append (ma.GetNEdges());
    else
      ndlevel.Append (ma.GetNFaces());
  }

  
  int RaviartThomasFESpace :: GetNDof () const
  {
    return ndlevel.Last();
  }
  
  int RaviartThomasFESpace :: GetNDofLevel (int level) const
  {
    return ndlevel[level];
  }
  
  
  
  void RaviartThomasFESpace :: GetDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    ARRAY<int> forient(6);
    
    if (ma.GetDimension() == 2)
      GetMeshAccess().GetElEdges (elnr, dnums, forient);
    else
      GetMeshAccess().GetElFaces (elnr, dnums, forient);
    
    
    if (!DefinedOn (ma.GetElIndex (elnr)))
      dnums = -1;

    
    // (*testout) << "el = " << elnr << ", dofs = " << dnums << endl;
  }
  
  
  void RaviartThomasFESpace :: GetSDofNrs (int selnr, ARRAY<int> & dnums) const
  {
    /*
      int eoa[12];
      ARRAY<int> eorient(12, eoa);
      GetMeshAccess().GetSElEdges (selnr, dnums, eorient);
      
      if (!DefinedOnBoundary (ma.GetSElIndex (selnr)))
      dnums = -1;
    */
    dnums.SetSize (1);
    dnums = -1;
    
    // (*testout) << "sel = " << selnr << ", dofs = " << dnums << endl;
  }

  Table<int> * RaviartThomasFESpace :: CreateSmoothingBlocks (int type) const
  {
    return 0;
  }
  
  
  void RaviartThomasFESpace :: 
  VTransformMR (int elnr, bool boundary,
		FlatMatrix<double> & mat, TRANSFORM_TYPE tt) const
  {
    int nd = 3;
    
    if (boundary) return;

    Vector<> fac(nd);
    
    GetTransformationFactors (elnr, fac);
    
    int i, j, k, l;
    
    if (tt & TRANSFORM_MAT_LEFT)
      for (k = 0; k < dimension; k++)
	for (i = 0; i < nd; i++)
	  for (j = 0; j < mat.Width(); j++)
	    mat(k+i*dimension, j) *= fac(i);
	
    if (tt & TRANSFORM_MAT_RIGHT)
      for (l = 0; l < dimension; l++)
	for (i = 0; i < mat.Height(); i++)
	  for (j = 0; j < nd; j++)
	    mat(i, l+j*dimension) *= fac(j);
  }
  
    
  void  RaviartThomasFESpace ::
  VTransformVR (int elnr, bool boundary,
		FlatVector<double> & vec, TRANSFORM_TYPE tt) const
  {
    int nd = 3;
    
    if (boundary) return;

    Vector<> fac(nd);
    
    GetTransformationFactors (elnr, fac);
    
    if ((tt & TRANSFORM_RHS) || (tt & TRANSFORM_SOL))
      {
	for (int k = 0; k < dimension; k++)
	  for (int i = 0; i < nd; i++)
	    vec(k+i*dimension) *= fac(i);
      }
  }
  
  void RaviartThomasFESpace ::
  GetTransformationFactors (int elnr, FlatVector<> & fac) const
  {
    ARRAY<int> edge_nums, edge_orient;
    
    fac = 1;

    ma.GetElEdges (elnr, edge_nums, edge_orient);
    for (int i = 0; i < 3; i++)
      fac(i) = edge_orient[i];
  }



  
  // register FESpaces
  namespace {
    class Init
    { 
    public: 
      Init ();
    };
    
    Init::Init()
    {
      GetFESpaceClasses().AddFESpace ("hdiv", RaviartThomasFESpace::Create);
    }
    
    Init init;
  }
}