/*********************************************************************/
/* File:   l2hofespace.cpp                                           */
/* Author: Start                                                     */
/* Date:   24. Feb. 2003                                             */
/*********************************************************************/

/**
   High Order Finite Element Space for L2
*/



#include <comp.hpp>
#include <multigrid.hpp>
using namespace ngmg;

namespace ngcomp
{
  using namespace ngcomp;

  L2HighOrderFESpace ::  
  L2HighOrderFESpace (const MeshAccess & ama, const Flags & flags)
    : FESpace (ama, 
	       int (flags.GetNumFlag ("order", 1)),
	       int (flags.GetNumFlag ("dim", 1)),
	       flags.GetDefineFlag ("complex"))
  {
    segm = new H1HighOrderSegm<IntegratedLegendreMonomialExt> (order);
    trig = new H1HighOrderTrig<IntegratedLegendreMonomialExt> (order);
    quad = new H1HighOrderQuad<IntegratedLegendreMonomialExt> (order);
    tet =  new H1HighOrderTet<IntegratedLegendreMonomialExt> (order);
    prism = new H1HighOrderPrism<IntegratedLegendreMonomialExt> (order);
    pyramid = new H1HighOrderPyramid<IntegratedLegendreMonomialExt> (order);
    hex = new H1HighOrderHex<IntegratedLegendreMonomialExt> (order);

    /*
      segm = new H1HighOrderSegm (order);
      trig = new H1HighOrderTrig<TrigExtensionMonomial> (order);
      quad = new H1HighOrderQuad<TrigExtensionMonomial> (order);
      tet =  new H1HighOrderTet<TrigExtensionMonomial> (order);
      prism = new H1HighOrderPrism<TrigExtensionMonomial> (order);
      pyramid = new H1HighOrderPyramid<TrigExtensionMonomial> (order);
      hex = new H1HighOrderHex<TrigExtensionMonomial> (order);
    */
    

    if (ma.GetDimension() == 2)
      {
	evaluator = 
	  new MassIntegrator<2> (new ConstantCoefficientFunction(1));
      }
    else
      {
	evaluator = 
	  new MassIntegrator<3> (new ConstantCoefficientFunction(1));
      }

    if (dimension > 1)
      {
	evaluator = new BlockBilinearFormIntegrator (*evaluator, dimension);
      }
  }


  L2HighOrderFESpace :: ~L2HighOrderFESpace ()
  {
    ;
  }


  FESpace * L2HighOrderFESpace :: 
  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 ElementFESpace (ma, order, dim, iscomplex);
    else
      return new L2HighOrderFESpace (ma, flags);
  }  

  void L2HighOrderFESpace :: Update()
  {
    nel = ma.GetNE();
    ndof = 0;
    first_element_dof.SetSize(nel+1);
    for (int i = 0; i < nel; i++)
      {
	first_element_dof[i] = ndof;
	switch (ma.GetElType(i))
	  {
	  case ET_TRIG:
	    ndof += (order+1)*(order+2)/2;
	    break;
	  case ET_QUAD:
	    ndof += (order+1)*(order+1);
	    break;
	  case ET_TET:
	    ndof += (order+1)*(order+2)*(order+3)/6;
	    break;
	  case ET_PRISM:
	    ndof += (order+1)*(order+2)*(order+1)/2;
	    break;
	  case ET_PYRAMID:
	    ndof += 5;
	    break;
	  case ET_HEX:
	    ndof += (order+1)*(order+1)*(order+1);
	    break;
	  }
      }
    first_element_dof[nel] = ndof;
  }



  const FiniteElement & L2HighOrderFESpace :: GetFE (int elnr, LocalHeap & lh) const
  {
    FiniteElement * fe = 0;
    
    switch (ma.GetElType(elnr))
      {
      case ET_TET:
	fe = tet; break;
      case ET_PYRAMID:
	fe = pyramid; break;
      case ET_PRISM:
	fe = prism; break;
      case ET_HEX:
	fe = hex; break;
      case ET_TRIG:
	fe = trig; break;
      case ET_QUAD:
	fe = quad; break;
      default:
	fe = 0;
      }
    
    ArrayMem<int,12> vnums; // calls GetElPNums -> max 12 for PRISM12
    ma.GetElVertices(elnr, vnums);
          
    if (!fe)
      {
	stringstream str;
	str << "L2HighOrderFESpace " << GetClassName() 
	    << ", undefined eltype " 
	    << ElementTopology::GetElementName(ma.GetElType(elnr))
	    << ", order = " << order << endl;
	throw Exception (str.str());
      }

    dynamic_cast<H1HighOrderFiniteElement*> (fe)
      -> SetVertexNumbers (vnums);

    return *fe;
  }
 
  const FiniteElement & L2HighOrderFESpace :: GetSFE (int elnr, LocalHeap & lh) const
  {
    FiniteElement * fe = 0;
    
    switch (ma.GetSElType(elnr))
      {
      case ET_SEGM:
	fe = segm; break;
      case ET_TRIG:
	fe = trig; break;
      case ET_QUAD:
	fe = quad; break;
      default:
	fe = 0;
      }
    
    ARRAY<int> vnums;
    ma.GetSElVertices(elnr, vnums);

    if (!fe)
      {
	stringstream str;
	str << "FESpace " << GetClassName() 
	    << ", undefined surface eltype " << ma.GetSElType(elnr) 
	    << ", order = " << order << endl;
	throw Exception (str.str());
      }

    dynamic_cast<H1HighOrderFiniteElement*> (fe)
      -> SetVertexNumbers (vnums);

    return *fe;
  }
 



  int L2HighOrderFESpace :: GetNDof () const
  {
    return ndof;
  }

  void L2HighOrderFESpace :: GetDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    int i, j;

    dnums.SetSize(0);
    int first = first_element_dof[elnr];
    int neldofs = first_element_dof[elnr+1] - first;
    for (j = 0; j < neldofs; j++)
      dnums.Append (first+j);

    if (!DefinedOn (ma.GetElIndex (elnr)))
      dnums = -1;
  }
  
  void L2HighOrderFESpace :: 
  GetSDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    dnums.SetSize (0);
  }
  
  Table<int> * L2HighOrderFESpace :: 
  CreateSmoothingBlocks ( int type) const
  {
    int i, j, first;
    ARRAY<int> cnt(nel);
    cnt = 0;
    for (i = 0; i < nel; i++)
      cnt[i] = first_element_dof[i+1]-first_element_dof[i];
	
    Table<int> & table = *new Table<int> (cnt);
    
    for (i = 0; i < nel; i++)
      {
	first = first_element_dof[i];
	for (j = 0; j < cnt[i]; j++)
	  table[i][j] = first+j;
      }
    return &table;
  }




















  L2SurfaceHighOrderFESpace ::  
  L2SurfaceHighOrderFESpace (const MeshAccess & ama, const Flags & flags)
    : FESpace (ama, 
	       int (flags.GetNumFlag ("order", 1)),
	       int (flags.GetNumFlag ("dim", 1)),
	       flags.GetDefineFlag ("complex"))
  {
    segm = new H1HighOrderSegm<IntegratedLegendreMonomialExt> (order);
    trig = new H1HighOrderTrig<IntegratedLegendreMonomialExt> (order);
    quad = new H1HighOrderQuad<IntegratedLegendreMonomialExt> (order);

    if (ma.GetDimension() == 2)
      {
	boundary_evaluator = 
	  new RobinIntegrator<2> (new ConstantCoefficientFunction(1));
      }
    else
      {
	boundary_evaluator = 
	  new RobinIntegrator<3> (new ConstantCoefficientFunction(1));
      }

    if (dimension > 1)
      {
	boundary_evaluator = 
	  new BlockBilinearFormIntegrator (*boundary_evaluator, dimension);
      }
  }


  L2SurfaceHighOrderFESpace :: ~L2SurfaceHighOrderFESpace ()
  {
    ;
  }

  FESpace * L2SurfaceHighOrderFESpace :: 
  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");

    return new L2SurfaceHighOrderFESpace (ma, flags);
  }  

  

  void L2SurfaceHighOrderFESpace :: Update()
  {
    int i;

    nel = ma.GetNSE();

    ndof = 0;
    first_element_dof.SetSize(nel+1);
    for (i = 0; i < nel; i++)
      {
	first_element_dof[i] = ndof;
	switch (ma.GetSElType(i))
	  {
	  case ET_SEGM:
	    ndof += order+1;
	    break;
	  case ET_TRIG:
	    ndof += (order+1)*(order+2)/2;
	    break;
	  case ET_QUAD:
	    ndof += (order+1)*(order+1);
	    break;
	  default:
	    ;
	  }
      }
    first_element_dof[nel] = ndof;
  }



  const FiniteElement & L2SurfaceHighOrderFESpace :: GetSFE (int elnr, LocalHeap & lh) const
  {
    FiniteElement * fe = 0;
    
    switch (ma.GetSElType(elnr))
      {
      case ET_SEGM:
	fe = segm; break;
      case ET_TRIG:
	fe = trig; break;
      case ET_QUAD:
	fe = quad; break;
      default:
	fe = 0;
      }
    
    ArrayMem<int,12> vnums; // calls GetElPNums -> max 12 for PRISM12
    ma.GetSElVertices(elnr, vnums);

    if (!fe)
      {
	stringstream str;
	str << "L2SurfaceHighOrderFESpace " << GetClassName() 
	    << ", undefined eltype " 
	    << ElementTopology::GetElementName(ma.GetSElType(elnr))
	    << ", order = " << order << endl;
	throw Exception (str.str());
      }

    dynamic_cast<H1HighOrderFiniteElement*> (fe)
      -> SetVertexNumbers (vnums);

    return *fe;
  }
 



  const FiniteElement & L2SurfaceHighOrderFESpace :: GetFE (int elnr, LocalHeap & lh) const
  {
    throw Exception ("Volume elements not available for L2SurfaceHighOrderFESpace");
  }
 


  int L2SurfaceHighOrderFESpace :: GetNDof () const
  {
    return ndof;
  }

  void L2SurfaceHighOrderFESpace :: GetSDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    int i, j;

    dnums.SetSize(0);
    int first = first_element_dof[elnr];
    int neldofs = first_element_dof[elnr+1] - first;
    for (j = 0; j < neldofs; j++)
      dnums.Append (first+j);

    if (!DefinedOn (ma.GetSElIndex (elnr)))
      dnums = -1;
    // (*testout) << "dnums = " << dnums << endl;
  }
  
  void L2SurfaceHighOrderFESpace :: 
  GetDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    dnums.SetSize (0);
  }
  
  Table<int> * L2SurfaceHighOrderFESpace :: 
  CreateSmoothingBlocks ( int type) const
  {
    int i, j, first;
    ARRAY<int> cnt(nel);
    cnt = 0;
    for (i = 0; i < nel; i++)
      cnt[i] = first_element_dof[i+1]-first_element_dof[i];
	
    Table<int> & table = *new Table<int> (cnt);
    
    for (i = 0; i < nel; i++)
      {
	first = first_element_dof[i];
	for (j = 0; j < cnt[i]; j++)
	  table[i][j] = first+j;
      }
    return &table;
  }



    
  // register FESpaces
  namespace {
  
    class Init
    { 
    public: 
      Init ();
    };
    
    Init::Init()
    {
      GetFESpaceClasses().AddFESpace ("l2", L2HighOrderFESpace::Create);
      GetFESpaceClasses().AddFESpace ("l2surf", L2SurfaceHighOrderFESpace::Create);
    }
    
    Init init;
  }


}