/*********************************************************************/
/* File:   h1hofespace.cpp                                           */
/* Author: Start                                                     */
/* Date:   10. Feb. 2003                                             */
/*********************************************************************/

/**
   High Order Finite Element Space
*/

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

namespace ngcomp
{
  using namespace ngcomp;


  H1HighOrderFESpace ::  
  H1HighOrderFESpace (const MeshAccess & ama, const Flags & flags)
    : FESpace (ama, flags)
  {
    //    order_inner = int (flags.GetNumFlag ("orderinner", order));

    augmented = int (flags.GetNumFlag ("augmented", 0));
    rel_order = int (flags.GetNumFlag ("relorder", 0));
    //tet_order = int (flags.GetNumFlag ("increase_tet_order",0)); 


    
    Flags loflags;
    loflags.SetFlag ("order", 1);
    loflags.SetFlag ("dim", dimension);
    if (iscomplex) loflags.SetFlag ("comples");
    low_order_space = new NodalFESpace (ma, loflags);

    // low_order_space = new NodalFESpace (ma, 1, dimension, iscomplex);
  
    uniform_order_inner = int (flags.GetNumFlag ("orderinner", order));

    /*
    if (flags.GetDefineFlag ("optext"))
      {
	segm = new H1HighOrderSegm<TrigExtensionOptimal> (order);
	trig = new H1HighOrderTrig<TrigExtensionOptimal> (order);
	tet =  new H1HighOrderTet<TrigExtensionOptimal> (order);
	prism = new H1HighOrderPrism<TrigExtensionOptimal> (order);
	pyramid = new H1HighOrderPyramid<TrigExtensionOptimal> (order);
	quad = new H1HighOrderQuad<TrigExtensionOptimal> (order);
	hex = new H1HighOrderHex<TrigExtensionOptimal> (order);
      }

    else if (flags.GetDefineFlag ("minext"))
      {
	segm = new H1HighOrderSegm<TrigExtensionMin> (order);
	trig = new H1HighOrderTrig<TrigExtensionMin> (order);
	tet =  new H1HighOrderTet<TrigExtensionMin> (order);
	prism = new H1HighOrderPrism<TrigExtensionMin> (order);
	pyramid = new H1HighOrderPyramid<TrigExtensionMin> (order);
	quad = new H1HighOrderQuad<TrigExtensionMin> (order);
	hex = new H1HighOrderHex<TrigExtensionMin> (order);
      }
    else
    */
      {
	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);
      }

    
    dynamic_cast<H1HighOrderFiniteElement*> (segm) -> SetAugmented (augmented);
    dynamic_cast<H1HighOrderFiniteElement*> (trig) -> SetAugmented (augmented);
    dynamic_cast<H1HighOrderFiniteElement*> (tet) -> SetAugmented (augmented);
    dynamic_cast<H1HighOrderFiniteElement*> (prism) -> SetAugmented (augmented);
    dynamic_cast<H1HighOrderFiniteElement*> (quad) -> SetAugmented (augmented);
    dynamic_cast<H1HighOrderFiniteElement*> (pyramid) -> SetAugmented (augmented);
    dynamic_cast<H1HighOrderFiniteElement*> (hex) -> SetAugmented (augmented);


    static ConstantCoefficientFunction one(1);
    if (ma.GetDimension() == 2)
      {
	evaluator = new MassIntegrator<2> (&one);
	boundary_evaluator = 0;
      }
    else
      {
	evaluator = new MassIntegrator<3> (&one);
	boundary_evaluator = new RobinIntegrator<3> (&one);
      }

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

    prol = new LinearProlongation(*this);
  }


  H1HighOrderFESpace :: ~H1HighOrderFESpace ()
  {
    ;
  }


  FESpace * H1HighOrderFESpace :: 
  Create (const MeshAccess & ma, const Flags & flags)
  {
    return new H1HighOrderFESpace (ma, flags);
  }
    

  void H1HighOrderFESpace :: Update()
  {
    // cout << "update h1ho space" << endl;
    low_order_space -> Update();

    nv = ma.GetNV();
    ned = ma.GetNEdges();
    nfa = ma.GetNFaces();
    nel = ma.GetNE();


    order_edge.SetSize (ned);
    order_face.SetSize (nfa);
    order_inner.SetSize (nel);
   
    /*
      order_edge = order;
      order_face = order;
      order_inner = order;
    */
    
      order_edge = 1;
      order_face = 1;
      order_inner = uniform_order_inner;
    
    

    ARRAY<int> eledges, elfaces, vnums;

    for (int i = 0; i < nel; i++)
      {
	if (!DefinedOn (ma.GetElIndex(i))) continue;

	int elorder = ma.GetElOrder(i) + rel_order;

	if (order > elorder) elorder = order;

	ma.GetElEdges (i, eledges);
	for (int j = 0; j < eledges.Size(); j++)
	  if (elorder > order_edge[eledges[j]])
	    order_edge[eledges[j]] = elorder;
	
	if (ma.GetDimension() == 3)
	  {
	    ma.GetElFaces (i, elfaces);
	    for (int j = 0; j < elfaces.Size(); j++)
	      {
		if (elorder > order_face[elfaces[j]])
		  order_face[elfaces[j]] = elorder;
		//if(ma.GetElType(i)==ET_TET) order_face[elfaces[j]] *=2; // TrigFaces Test SZ 
	      }
	  }

	  
	if (elorder > order_inner[i])
	  order_inner[i] = elorder;
	//if(ma.GetElType(i)==ET_TET) order_inner[i] *=2; // TetInner *2 Test SZ  
      }


    // order_inner = 1;

    /*
    // skip dofs on Dirichlet boundary hack [JS]
    int nsel = ma.GetNE();
    for (i = 0; i < nsel; i++)
      {
	  ma.GetSElEdges (i, eledges);
	  int face = ma.GetSElFace (i);
	  
	  for (j = 0; j < eledges.Size(); j++)
	    order_edge[eledges[j]] = 1;
	  order_face[face] = 1;
      }
    */



    /*
    (*testout) << "order_edge = " << order_edge << endl;
    if(ma.GetDimension() == 3) (*testout) << "order_face = " << order_face << endl;
    (*testout) << "order_inner = " << order_inner << endl;
    */
    UpdateDofTables ();
  }


  void H1HighOrderFESpace :: UpdateDofTables ()
  {
    nv = ma.GetNV();
    ned = ma.GetNEdges();
    nfa = ma.GetNFaces();
    nel = ma.GetNE();

    ndof = nv;
    if (augmented == 1) ndof += nv;
    if (augmented == 2) ndof = order*nv;

    first_edge_dof.SetSize (ned+1);
    for (int i = 0; i < ned; i++)
      {
	first_edge_dof[i] = ndof;
	ndof += order_edge[i] - 1;
      }
    first_edge_dof[ned] = ndof;
    
    first_face_dof.SetSize (nfa+1);
    ARRAY<int> fapnums;
    for (int i = 0; i < nfa; i++)
      {
	ma.GetFacePNums (i, fapnums);
	first_face_dof[i] = ndof;
	int p = order_face[i];
	if (fapnums.Size() == 3)
	  ndof += (p-1)*(p-2)/2;
	else
	  ndof += (p-1)*(p-1);
      }
    first_face_dof[nfa] = ndof;
 
    first_element_dof.SetSize(nel+1);
    for (int i = 0; i < nel; i++)
      {
	first_element_dof[i] = ndof;
	int p = order_inner[i];	
	switch (ma.GetElType(i))
	  {
	  case ET_TRIG:
	    ndof += (p-1)*(p-2)/2;
	    break;
	  case ET_QUAD:
	    ndof += (p-1)*(p-1);
	    break;
	  case ET_TET:
	    ndof += (p-1)*(p-2)*(p-3)/6;
	    break;
	  case ET_PRISM:
	    ndof += (p-1)*(p-2)*(p-1)/2;
	    break;
	  case ET_PYRAMID:
	    ndof += (p-1)*(p-2)*(2*p-3)/6;
	    break;
	  case ET_HEX:
	    ndof += (p-1)*(p-1)*(p-1);
	    break;
	  }
      }
    first_element_dof[nel] = ndof;
    
    /*
    (*testout) << "augmented " << augmented << endl; 
    (*testout) << "first edge = " << first_edge_dof << endl;
    (*testout) << "first face = " << first_face_dof << endl;
    (*testout) << "first inner = " << first_element_dof << endl;
    */

    while (ma.GetNLevels() > ndlevel.Size())
      ndlevel.Append (ndof);
    ndlevel.Last() = ndof;
    prol->Update();
  }


  void H1HighOrderFESpace :: PrintReport (ostream & ost)
  {
    FESpace::PrintReport (ost);

//     (*testout) << "first edge = " << first_edge_dof << endl;
//     (*testout) << "first face = " << first_face_dof << endl;
//     (*testout) << "first inner = " << first_element_dof << endl;

  }


  const FiniteElement & H1HighOrderFESpace :: GetFE (int elnr, LocalHeap & lh) const
  {
	  try
	  {
    /*
    // this version is not thead-safe, but elements with different orthogonal-pols

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

    H1HighOrderFiniteElement * hofe =
      dynamic_cast<H1HighOrderFiniteElement*> (fe);
    */


    H1HighOrderFiniteElement * hofe;
    switch (ma.GetElType(elnr))
      {
      case ET_TET:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderTet<>)))  H1HighOrderTet<> (order);
	  break;
	}
      case ET_PYRAMID:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderPyramid<>)))  H1HighOrderPyramid<> (order);
	  break;
	}
      case ET_PRISM:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderPrism<>)))  H1HighOrderPrism<> (order);
	  break;
	}
      case ET_HEX:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderHex<>)))  H1HighOrderHex<> (order);
	  break;
	}
      case ET_TRIG:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderTrig<>)))  H1HighOrderTrig<> (order);
	  break;
	}
      case ET_QUAD:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderQuad<>)))  H1HighOrderQuad<> (order);
	  break;
	}
      default:
	hofe = 0; 
      }


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



    hofe -> SetAugmented (augmented);
    
    ArrayMem<int,12> vnums; // calls GetElPNums -> max 12 for PRISM12
    ArrayMem<int, 12> ednums, order_ed;
    ArrayMem<int, 6> fanums, order_fa;

    ma.GetElVertices(elnr, vnums);
    ma.GetElEdges(elnr, ednums);

    order_ed.SetSize (ednums.Size());

    for (int j = 0; j < ednums.Size(); j++)
      order_ed[j] = order_edge[ednums[j]];

    hofe -> SetVertexNumbers (vnums);
    hofe -> SetOrderEdge (order_ed);
    
    if (ma.GetDimension() == 3)
      {
	ma.GetElFaces(elnr, fanums);
	order_fa.SetSize (fanums.Size());
	for (int j = 0; j < fanums.Size(); j++)
	  order_fa[j] = order_face[fanums[j]];
	
	hofe -> SetOrderFace (order_fa);
      }


    hofe -> SetOrderInner (order_inner[elnr]);

    return *hofe;
	}
	    catch (Exception & e)
      {
		  e.Append ("in H1HoFESpace::GetElement");
	e.Append ("\n");
	throw;
      }
  }
 



  const FiniteElement & H1HighOrderFESpace :: 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;
      }

    H1HighOrderFiniteElement * hofe =
      dynamic_cast<H1HighOrderFiniteElement*> (fe);
    */


    H1HighOrderFiniteElement * hofe;
    switch (ma.GetSElType(elnr))
      {
      case ET_TRIG:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderTrig<>)))  H1HighOrderTrig<> (1);
	  break;
	}
      case ET_QUAD:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderQuad<>)))  H1HighOrderQuad<> (1);
	  break;
	}
      case ET_SEGM:
	{
	  hofe = new (lh.Alloc (sizeof(H1HighOrderSegm<>)))  H1HighOrderSegm<> (1);
	  break;
	}
      default:
	hofe = 0; 
      }

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

    hofe -> SetAugmented (augmented);
    
    ArrayMem<int, 8> vnums;
    ArrayMem<int, 4> ednums, order_ed;

    ma.GetSElVertices(elnr, vnums);
    ma.GetSElEdges(elnr, ednums);

    order_ed.SetSize (ednums.Size());

    for (int j = 0; j < ednums.Size(); j++)
      order_ed[j] = order_edge[ednums[j]];


    hofe -> SetVertexNumbers (vnums);

    if (ma.GetDimension() == 3)
      {
	hofe -> SetOrderEdge (order_ed);
	hofe -> SetOrderInner (order_face[ma.GetSElFace(elnr)]);
      }
    else
      {
	hofe -> SetOrderInner (order_ed[0]);
      }
    return *hofe;
  }
 






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



  int H1HighOrderFESpace :: GetNDofLevel (int level) const
  {
    return ndlevel[level];
  }



  void H1HighOrderFESpace :: GetDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    ArrayMem<int,12> vnums, ednums, fanums; // , eorient, forient;
    int i, j;

    ma.GetElVertices (elnr, vnums);
    ma.GetElEdges (elnr, ednums);
    if (ma.GetDimension() == 3)
      ma.GetElFaces (elnr, fanums);


    // (*testout) << "vnums = " << vnums << endl
    //    (*testout) << "ednums = " << ednums << endl;
      // << "fanums = " << fanums << endl;

    dnums.SetSize(0);
    for (i = 0; i < vnums.Size(); i++)
      dnums.Append (vnums[i]);

    if (augmented == 1)
      for (i = 0; i < vnums.Size(); i++)
	dnums.Append (nv+vnums[i]);
    else if (augmented == 2)
      for (i = 0; i < vnums.Size(); i++)
	for (j = 0; j < order-1; j++)
	  dnums.Append (nv+(order-1)*vnums[i]+j);

    for (i = 0; i < ednums.Size(); i++)
      {
	int first = first_edge_dof[ednums[i]];
	int neddofs = first_edge_dof[ednums[i]+1] - first;
	// (*testout) << "first = " << first << "; ned = " << neddofs << endl;
	for (j = 0; j < neddofs; j++)
	  dnums.Append (first+j);
      }

    if (ma.GetDimension() == 3)
      for (i = 0; i < fanums.Size(); i++)
	{
	  int first = first_face_dof[fanums[i]];
	  int nfadofs = first_face_dof[fanums[i]+1] - first;
	  // (*testout) << "first = " << first << "; nfa = " << nfadofs << endl;
	  for (j = 0; j < nfadofs; j++)
	    dnums.Append (first+j);
	}

    // (*testout) << "has faces" << endl;


    int first = first_element_dof[elnr];
    int neldofs = first_element_dof[elnr+1] - first;

    //    if (!eliminate_internal)
      for (j = 0; j < neldofs; j++)
	dnums.Append (first+j);
      /*
    else
      for (j = 0; j < neldofs; j++)
	dnums.Append (-1);
      */

    if (!DefinedOn (ma.GetElIndex (elnr)))
      dnums = -1;

    /*
    (*testout) << "getdnums, el = " << elnr << endl 
	       << "vnums = " << vnums << endl
	       << "ednums = " << ednums << ", fanums = " << fanums << endl;
    // (*testout) << "dnums(" << elnr << ") = " << endl << dnums << endl;
    */
  }



  void  H1HighOrderFESpace :: GetExternalDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    if (!eliminate_internal) 
      {
	GetDofNrs (elnr, dnums);
	return;
      }

    ArrayMem<int,12> vnums, ednums, fanums; 
    int i, j;

    ma.GetElVertices (elnr, vnums);
    ma.GetElEdges (elnr, ednums);
    if (ma.GetDimension() == 3)
      ma.GetElFaces (elnr, fanums);

    dnums.SetSize(0);
    for (i = 0; i < vnums.Size(); i++)
      dnums.Append (vnums[i]);

    if (augmented == 1)
      for (i = 0; i < vnums.Size(); i++)
	dnums.Append (nv+vnums[i]);
    else if (augmented == 2)
      for (i = 0; i < vnums.Size(); i++)
	for (j = 0; j < order-1; j++)
	  dnums.Append (nv+(order-1)*vnums[i]+j);

    for (i = 0; i < ednums.Size(); i++)
      {
	int first = first_edge_dof[ednums[i]];
	int neddofs = first_edge_dof[ednums[i]+1] - first;
	for (j = 0; j < neddofs; j++)
	  dnums.Append (first+j);
      }

    if (ma.GetDimension() == 3)
      for (i = 0; i < fanums.Size(); i++)
	{
	  int first = first_face_dof[fanums[i]];
	  int nfadofs = first_face_dof[fanums[i]+1] - first;
	  for (j = 0; j < nfadofs; j++)
	    dnums.Append (first+j);
	}

    if (!DefinedOn (ma.GetElIndex (elnr)))
      dnums = -1;
  }



  
  void H1HighOrderFESpace :: 
  GetSDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    ARRAY<int> vnums, ednums;
    int fanum; // , faorient;
    int i, j;

    ma.GetSElVertices (elnr, vnums);
    ma.GetSElEdges (elnr, ednums);
    if (ma.GetDimension() == 3)
      fanum = ma.GetSElFace (elnr); // , fanum, faorient);

    dnums.SetSize(0);
    for (i = 0; i < vnums.Size(); i++)
      dnums.Append (vnums[i]);

    if (augmented == 1)
      for (i = 0; i < vnums.Size(); i++)
	dnums.Append (nv+vnums[i]);
    else if (augmented == 2)
      for (i = 0; i < vnums.Size(); i++)
	for (j = 0; j < order-1; j++)
	  dnums.Append (nv+(order-1)*vnums[i]+j);

    for (i = 0; i < ednums.Size(); i++)
      {
	int first = first_edge_dof[ednums[i]];
	int neddofs = first_edge_dof[ednums[i]+1] - first;
	for (j = 0; j < neddofs; j++)
	  dnums.Append (first+j);
	// dnums.Append (-1);
      }

    if (ma.GetDimension() == 3)
      {
	int first = first_face_dof[fanum];
	int nfadofs = first_face_dof[fanum+1] - first;
	for (j = 0; j < nfadofs; j++)
	  //	  dnums.Append (-1);
	  dnums.Append (first+j);
      }

    if (!DefinedOnBoundary (ma.GetSElIndex (elnr)))
      dnums = -1;

    //    (*testout) << "bound-dnums(" << elnr << ") = " << endl << dnums << endl;
  }
  
  Table<int> * H1HighOrderFESpace :: 
  CreateSmoothingBlocks ( int type) const
  {
    int i, j, k;

    if (0)
      {
	ARRAY<int> cnt(2);
	cnt = 0;
	cnt[0] = 2*nv;
	cnt[1] = GetNDof()-2*nv;

	Table<int> & table = *new Table<int> (cnt); 
	
	cnt = 0;
	
	for (i = 0; i < 2*nv; i++)
	  {
	    table[0][cnt[0]++] = i;
	  }
	for ( ; i < GetNDof(); i++)
	  table[1][cnt[1]++] = i;

	// (*testout) << "table = " << endl << table << endl;
	return &table;
	/*
	ARRAY<int> cnt(4);
	cnt = 0;
	cnt[0] = nv;
	for (i = 0; i < ned; i++)
	  cnt[1] += first_edge_dof[i+1] - first_edge_dof[i];
	for (i = 0; i < nfa; i++)
	  cnt[1] += first_face_dof[i+1] - first_face_dof[i];
	for (i = 0; i < nel; i++)
	  cnt[1] += first_element_dof[i+1] - first_element_dof[i];


	Table<int> & table = *new Table<int> (cnt); 
	
	cnt = 0;
	
	for (i = 0; i < nv; i++)
	  {
	    table[0][cnt[0]++] = i;
	  }
	for (i = 0; i < ned; i++)
	  {
	    for (j = first_edge_dof[i]; j < first_edge_dof[i+1]; j++)
	      table[1][cnt[1]++] = j;
	  }
	for (i = 0; i < nfa; i++)
	  {
	    for (j = first_face_dof[i]; j < first_face_dof[i+1]; j++)
	      table[1][cnt[1]++] = j;      
	  }
	for (i = 0; i < nel; i++)
	  for (j = first_element_dof[i]; j < first_element_dof[i+1]; j++)
	    table[1][cnt[1]++] = j;      

	(*testout) << "table = " << endl << table << endl;
	return &table;
	*/
      }


    if (ma.GetDimension() == 3)
      {
	
	  // V + E + F + I
	int i, j, first;
	ARRAY<int> cnt(nv+ned+nfa+nel);
	cnt = 0;

	int nvdof = 1;
	if (augmented == 1) nvdof = 2;
	if (augmented == 2) nvdof = order;
	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;
	for (i = 0; i < ned; i++)
	  cnt[nv+i] = first_edge_dof[i+1]-first_edge_dof[i];
	for (i = 0; i < nfa; i++)
	  cnt[nv+ned+i] = first_face_dof[i+1]-first_face_dof[i];
	if (!eliminate_internal)
	  for (i = 0; i < nel; i++)
	    cnt[nv+ned+nfa+i] = first_element_dof[i+1]-first_element_dof[i];
	else
	  for (i = 0; i < nel; i++)
	    cnt[nv+ned+nfa+i] = 0;
	  

	Table<int> & table = *new Table<int> (cnt);

	for (i = 0; i < nv; i++)
	  table[i][0] = i;
	if (augmented == 1)
	  for (i = 0; i < nv; i++)
	    table[i][1] = nv+i;
	else if (augmented == 2)
	  for (i = 0; i < nv; i++)
	    for (j = 0; j < order-1; j++)
	      table[i][j+1] = nv+i*(order-1)+j;
	for (i = 0; i < ned; i++)
	  {
	    first = first_edge_dof[i];
	    for (j = 0; j < cnt[nv+i]; j++)
	      table[nv+i][j] = first+j;
	  }
	for (i = 0; i < nfa; i++)
	  {
	    first = first_face_dof[i];
	    for (j = 0; j < cnt[nv+ned+i]; j++)
	      table[nv+ned+i][j] = first+j;
	  }
	if (!eliminate_internal)
	  for (i = 0; i < nel; i++)
	    {
	      first = first_element_dof[i];
	      for (j = 0; j < cnt[nv+ned+nfa+i]; j++)
		table[nv+ned+nfa+i][j] = first+j;
	    }
	// (*testout) << "H1HO-table = " << table << endl;
	return &table;

	


	/*
	// V + E + FI
	int i, j, first;
	ARRAY<int> cnt(nv+ned+nfa);
	ARRAY<int> vnums; 
	ARRAY<int> orient; 
	ARRAY<int> ednums, fanums, faorient;
	
	
	cnt = 0;
	int nvdof = 1;
	if (augmented == 1) nvdof = 2;
	if (augmented == 2) nvdof = order;

	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;
	for (i = 0; i < ned; i++)
	  cnt[nv+i] = first_edge_dof[i+1]-first_edge_dof[i];
	for (i = 0; i < nfa; i++)
	  cnt[nv+ned+i] = first_face_dof[i+1]-first_face_dof[i];
	for (i = 0; i < nel; i++)
	  {
	    ma.GetElFaces (i,fanums,orient);
	    for (k = 0; k < fanums.Size(); k++)
	      cnt[nv+ned+fanums[k]] += first_element_dof[i+1] - first_element_dof[i];
	  }

	Table<int> & table = *new Table<int> (cnt);
	cnt = 0;

	for (i = 0; i < nv; i++)
 	  table[i][0] = i;
	if (augmented == 1)
	  for (i = 0; i < nv; i++)
	    table[i][1] = nv+i;
	else if (augmented == 2)
	  for (i = 0; i < nv; i++)
	    for (j = 0; j < order-1; j++)
	      table[i][j+1] = nv+i*(order-1)+j;
	  
	for (i = 0; i < ned; i++)
	  {
	    first = first_edge_dof[i];
	    for (j = first; j < first_edge_dof[i+1]; j++)
	      table[nv+i][cnt[nv+i]++] = j;
	  }
	for (i = 0; i < nfa; i++)
	  {
	    first = first_face_dof[i];
	    for (j = first; j < first_face_dof[i+1]; j++)
	      table[nv+ned+i][cnt[nv+ned+i]++] = j;
	  }
	for (i = 0; i < nel; i++)
	  {
	    ma.GetElFaces (i,fanums,orient);
	    for (k = 0; k < fanums.Size(); k++)
	      {
		int dof = nv+ned+fanums[k];
		for (j = first_element_dof[i]; j < first_element_dof[i+1]; j++)
		  table[dof][cnt[dof]++] = j;
	      }
	  }
//	(*testout) << "H1HO-table = " << table << endl;
	return &table;
	*/

	/* 
	// V + EI + FI
	int i, j, first;
	ARRAY<int> cnt(nv+ned+nfa);
	ARRAY<int> vnums; 
	ARRAY<int> orient; 
	ARRAY<int> ednums, fanums, faorient;
	
	
	cnt = 0;
	int nvdof = 1;
	if (augmented == 1) nvdof = 2;
	if (augmented == 2) nvdof = order;

	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;
	for (i = 0; i < ned; i++)
	  cnt[nv+i] = first_edge_dof[i+1]-first_edge_dof[i];
	for (i = 0; i < nfa; i++)
	  cnt[nv+ned+i] = first_face_dof[i+1]-first_face_dof[i];
	for (i = 0; i < nel; i++)
	  {
	    ma.GetElEdges (i,ednums,orient);
	    for (k = 0; k < ednums.Size(); k++)
	      cnt[nv+ednums[k]] += first_element_dof[i+1] - first_element_dof[i];
	    ma.GetElFaces (i,fanums,orient);
	    for (k = 0; k < fanums.Size(); k++)
	      cnt[nv+ned+fanums[k]] += first_element_dof[i+1] - first_element_dof[i];
	  }

	Table<int> & table = *new Table<int> (cnt);
	cnt = 0;

	for (i = 0; i < nv; i++)
 	  table[i][0] = i;
	if (augmented == 1)
	  for (i = 0; i < nv; i++)
	    table[i][1] = nv+i;
	else if (augmented == 2)
	  for (i = 0; i < nv; i++)
	    for (j = 0; j < order-1; j++)
	      table[i][j+1] = nv+i*(order-1)+j;
	  
	for (i = 0; i < ned; i++)
	  {
	    first = first_edge_dof[i];
	    for (j = first; j < first_edge_dof[i+1]; j++)
	      table[nv+i][cnt[nv+i]++] = j;
	  }
	for (i = 0; i < nfa; i++)
	  {
	    first = first_face_dof[i];
	    for (j = first; j < first_face_dof[i+1]; j++)
	      table[nv+ned+i][cnt[nv+ned+i]++] = j;
	  }
	for (i = 0; i < nel; i++)
	  {
	    ma.GetElEdges (i,ednums,orient);
	    for (k = 0; k < ednums.Size(); k++)
	      {
		int dof = nv+ednums[k];
		for (j = first_element_dof[i]; j < first_element_dof[i+1]; j++)
		  table[dof][cnt[dof]++] = j;
	      }
	    ma.GetElFaces (i,fanums,orient);
	    for (k = 0; k < fanums.Size(); k++)
	      {
		int dof = nv+ned+fanums[k];
		for (j = first_element_dof[i]; j < first_element_dof[i+1]; j++)
		  table[dof][cnt[dof]++] = j;
	      }
	  }
	(*testout) << "H1HO-table = " << table << endl;
	return &table;
	
	*/ 

	/*

	
	// VE + F + I
	int i, j, k, first;
	ARRAY<int> cnt(nv+nfa+nel);
	cnt = 0;

	int nvdof = 1;
	if (augmented == 1) nvdof = 2;
	if (augmented == 2) nvdof = order;
	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;

	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    int ndof = first_edge_dof[i+1]-first_edge_dof[i];
	    ma.GetEdgePNums (i, v1, v2);
	    cnt[v1] += ndof;
	    cnt[v2] += ndof;
	  }

	for (i = 0; i < nfa; i++)
	  cnt[nv+i] = first_face_dof[i+1]-first_face_dof[i];

	if (!eliminate_internal)
	  for (i = 0; i < nel; i++)
	    cnt[nv+nfa+i] = first_element_dof[i+1]-first_element_dof[i];

	Table<int> & table = *new Table<int> (cnt);

	for (i = 0; i < nv; i++)
	  table[i][0] = i;
	if (augmented == 1)
	  for (i = 0; i < nv; i++)
	    table[i][1] = nv+i;
	else if (augmented == 2)
	  for (i = 0; i < nv; i++)
	    for (j = 0; j < order-1; j++)
	      table[i][j+1] = nv+i*(order-1)+j;
	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;

	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    first = first_edge_dof[i];
	    int ndof = first_edge_dof[i+1]-first;
	    ma.GetEdgePNums (i, v1, v2);
	    for (j = 0; j < ndof; j++)
	      {
		table[v1][cnt[v1]++] = first+j;
		table[v2][cnt[v2]++] = first+j;
	      }
	  }
	for (i = 0; i < nfa; i++)
	  {
	    first = first_face_dof[i];
	    int ndof = first_face_dof[i+1]-first_face_dof[i];
	    for (j = 0; j < ndof; j++)
	      table[nv+i][j] = first+j;
	    cnt[nv+i] = ndof;
	  }

	if (!eliminate_internal)
	  for (i = 0; i < nel; i++)
	    {
	      first = first_element_dof[i];
	      for (j = 0; j < cnt[nv+ned+nfa+i]; j++)
		table[nv+nfa+i][j] = first+j;
	    }

	    return &table;
	
	*/
	/*
	
	// V + EF + I
	int i, j, k, first;
	ARRAY<int> cnt(nv+ned+nel);
	ARRAY<int> f2ed;
	cnt = 0;

	int nvdof = 1;
	if (augmented == 1) nvdof = 2;
	if (augmented == 2) nvdof = order;
	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;

	for (i = 0; i < ned; i++)
	  cnt[nv+i]= first_edge_dof[i+1]-first_edge_dof[i];

	for (i = 0; i < nfa; i++)
	  {
	    int ndof = first_face_dof[i+1]-first_face_dof[i];
	    ma.GetFaceEdges (i, f2ed);
	    for (j = 0; j < f2ed.Size(); j++)
	      cnt[nv+f2ed[j]] += ndof;
	  }

	if (!eliminate_internal)
	  for (i = 0; i < nel; i++)
	    cnt[nv+ned+i] = first_element_dof[i+1]-first_element_dof[i];

	Table<int> & table = *new Table<int> (cnt);

	for (i = 0; i < nv; i++)
	  table[i][0] = i;
	if (augmented == 1)
	  for (i = 0; i < nv; i++)
	    table[i][1] = nv+i;
	else if (augmented == 2)
	  for (i = 0; i < nv; i++)
	    for (j = 0; j < order-1; j++)
	      table[i][j+1] = nv+i*(order-1)+j;

	cnt = 0;
	for (i = 0; i < ned; i++)
	  {
	    first = first_edge_dof[i];
	    int ndof = first_edge_dof[i+1]-first;
	    for (j = 0; j < ndof; j++)
	      table[nv+i][cnt[nv+i]++] = first+j;
	  }
	for (i = 0; i < nfa; i++)
	  {
	    first = first_face_dof[i];
	    int ndof = first_face_dof[i+1]-first;
	    ma.GetFaceEdges (i, f2ed);
	    for (k = 0; k < f2ed.Size(); k++)
	      for (j = 0; j < ndof; j++)
		table[nv+f2ed[k]][cnt[nv+f2ed[k]]++] = first+j;

	  }
	
	if (!eliminate_internal)
	  for (i = 0; i < nel; i++)
	    {
	      first = first_element_dof[i];
	      int ndof = first_element_dof[i+1]-first;
	      for (j = 0; j < ndof; j++)
		table[nv+ned+i][j] = first+j;
	    }

	(*testout) << "H1HO-table (V+EF) = " << endl << table << endl;
	return &table;
	
	*/
	
	/*
	// VE + FI
	int i, j, k, first;
	ARRAY<int> fanums, faorient;
	ARRAY<int> cnt(nv+ned+nfa+nel);
	cnt = 0;
	for (i = 0; i < nv; i++)
	  cnt[i] = 1;
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    int ndof = first_edge_dof[i+1]-first_edge_dof[i];
	    ma.GetEdgePNums (i, v1, v2);
	    cnt[v1] += ndof;
	    cnt[v2] += ndof;
	  }
	for (i = 0; i < nfa; i++)
	  cnt[nv+ned+i] = first_face_dof[i+1]-first_face_dof[i];

	for (i = 0; i < nel; i++)
	  {
	    ma.GetElFaces (i, fanums, faorient);
	    int ndof = first_element_dof[i+1] - first_element_dof[i];
	    for (j = 0; j < fanums.Size(); j++)
	      cnt[nv+ned+fanums[j]] += ndof;
	  }
	
	Table<int> & table = *new Table<int> (cnt);

	for (i = 0; i < nv; i++)
	  {
	    table[i][0] = i;
	    cnt[i] = 1;
	  }
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    first = first_edge_dof[i];
	    int ndof = first_edge_dof[i+1]-first;
	    ma.GetEdgePNums (i, v1, v2);
	    for (j = 0; j < ndof; j++)
	      {
		table[v1][cnt[v1]++] = first+j;
		table[v2][cnt[v2]++] = first+j;
	      }
	  }
	for (i = 0; i < nfa; i++)
	  {
	    first = first_face_dof[i];
	    int ndof = first_face_dof[i+1]-first_face_dof[i];
	    for (j = 0; j < ndof; j++)
	      table[nv+ned+i][j] = first+j;
	    cnt[nv+ned+i] = ndof;
	  }

	for (i = 0; i < nel; i++)
	  {
	    ma.GetElFaces (i, fanums, faorient);
	    first = first_element_dof[i];
	    int ndof = first_element_dof[i+1] - first_element_dof[i];
	    for (j = 0; j < fanums.Size(); j++)
	      for (k = 0; k < ndof; k++)
		table[nv+ned+fanums[j]][cnt[nv+ned+fanums[j]]++] = first + k;
	    // cnt[nv+ned+nfa+i] = first_element_dof[i+1]-first_element_dof[i];
	  }
	return &table;

	*/


	/*
	// VEF + I
	int i, j, k, first;
	ARRAY<int> fanums, faorient;
	ARRAY<int> cnt(nv+ned+nfa+nel);
	cnt = 0;
	for (i = 0; i < nv; i++)
	  cnt[i] = 1;
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    int ndof = first_edge_dof[i+1]-first_edge_dof[i];
	    ma.GetEdgePNums (i, v1, v2);
	    cnt[v1] += ndof;
	    cnt[v2] += ndof;
	  }
	for (i = 0; i < nfa; i++)
	  { 
	    ARRAY<int>  pnums; 
	    ma.GetFacePNums(i,pnums); 
	    int ndof =  first_face_dof[i+1] - first_face_dof[i];
	    for(j=0;j<pnums.Size();j++) 
	      cnt[pnums[j]] += ndof ;
	  }
	
	for (i = 0; i < nel; i++)
	  {
	    ARRAY<int>  pnums; 
	    int ndof = first_element_dof[i+1] - first_element_dof[i];
	    cnt[nv + ned + nfa + i ] += ndof;
	  }
	
	Table<int> & table = *new Table<int> (cnt);

	for (i = 0; i < nv; i++)
	  {
	    table[i][0] = i;
	    cnt[i] = 1;
	  }
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    first = first_edge_dof[i];
	    int ndof = first_edge_dof[i+1]-first;
	    ma.GetEdgePNums (i, v1, v2);
	    for (j = 0; j < ndof; j++)
	      {
		table[v1][cnt[v1]++] = first+j;
		table[v2][cnt[v2]++] = first+j;
	      }
	  }
	for (i = 0; i < nfa; i++)
	  {
	    ARRAY<int> pnums; 
	    ma.GetFacePNums(i,pnums); 
	    first = first_face_dof[i];
	    int ndof = first_face_dof[i+1]-first_face_dof[i];
	    for(k=0;k<pnums.Size();k++)
	      for (j = 0; j < ndof; j++)
		table[pnums[k]][cnt[pnums[k]]++] = first+j;
	  }

	for (i = 0; i < nel; i++)
	  {
	    first = first_element_dof[i];
	    int ndof = first_element_dof[i+1]-first;
	    for (j = 0; j < ndof; j++)
	      table[nv+ned+nfa+i][j] = first+j;
	  }
	
      
	return &table;
	*/
       
	/*

	// VEFI
	int i, j, k, first;
	ARRAY<int> fanums, faorient;
	ARRAY<int> cnt(nv+ned+nfa+nel);
	cnt = 0;
	for (i = 0; i < nv; i++)
	  cnt[i] = 1;
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    int ndof = first_edge_dof[i+1]-first_edge_dof[i];
	    ma.GetEdgePNums (i, v1, v2);
	    cnt[v1] += ndof;
	    cnt[v2] += ndof;
	  }
	for (i = 0; i < nfa; i++)
	  { 
	    ARRAY<int>  pnums; 
	    ma.GetFacePNums(i,pnums); 
	    int ndof =  first_face_dof[i+1] - first_face_dof[i];
	    for(j=0;j<pnums.Size();j++) 
	      cnt[pnums[j]] += ndof ;
	  }
	for (i = 0; i < nel; i++)
	  {
	    ARRAY<int>  pnums; 
	    ma.GetElPNums(i,pnums); 
	    int ndof = first_element_dof[i+1] - first_element_dof[i];
	    for (j = 0; j < pnums.Size(); j++)
	      cnt[pnums[j]] += ndof;
	  }
	
	Table<int> & table = *new Table<int> (cnt);

	for (i = 0; i < nv; i++)
	  {
	    table[i][0] = i;
	    cnt[i] = 1;
	  }
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    first = first_edge_dof[i];
	    int ndof = first_edge_dof[i+1]-first;
	    ma.GetEdgePNums (i, v1, v2);
	    for (j = 0; j < ndof; j++)
	      {
		table[v1][cnt[v1]++] = first+j;
		table[v2][cnt[v2]++] = first+j;
	      }
	  }
	for (i = 0; i < nfa; i++)
	  {
	    ARRAY<int> pnums; 
	    ma.GetFacePNums(i,pnums); 
	    first = first_face_dof[i];
	    int ndof = first_face_dof[i+1]-first_face_dof[i];
	    for(k=0;k<pnums.Size();k++)
	      for (j = 0; j < ndof; j++)
		table[pnums[k]][cnt[pnums[k]]++] = first+j;
	  }

	for (i = 0; i < nel; i++)
	  {
	    ARRAY<int> pnums; 
	    ma.GetElPNums(i,pnums); 
	    first = first_element_dof[i];
	    int ndof = first_element_dof[i+1]-first_element_dof[i];
	    for(k=0;k<pnums.Size();k++)
	      for (j = 0; j < ndof; j++)
		table[pnums[k]][cnt[pnums[k]]++] = first+j;
	  }
      
	return &table;
	
	*/

	
      }
    else
      {
	  // V + E + I
	int i, j, first;
	ARRAY<int> cnt(nv+ned+nel);
	cnt = 0;

	int nvdof = 1;
	if (augmented == 1) nvdof = 2;
	if (augmented == 2) nvdof = order;
	for (i = 0; i < nv; i++)
	  cnt[i] = nvdof;
	for (i = 0; i < ned; i++)
	  cnt[nv+i] = first_edge_dof[i+1]-first_edge_dof[i];
	for (i = 0; i < nel; i++)
	  cnt[nv+ned+i] = first_element_dof[i+1]-first_element_dof[i];
	
	Table<int> & table = *new Table<int> (cnt);
	
	for (i = 0; i < nv; i++)
	  table[i][0] = i;
	if (augmented == 1)
	  for (i = 0; i < nv; i++)
	    table[i][1] = nv+i;
	else if (augmented == 2)
	  for (i = 0; i < nv; i++)
	    for (j = 0; j < order-1; j++)
	      table[i][j+1] = nv+i*(order-1)+j;
	for (i = 0; i < ned; i++)
	  {
	    first = first_edge_dof[i];
	    for (j = 0; j < cnt[nv+i]; j++)
	      table[nv+i][j] = first+j;
	  }
	for (i = 0; i < nel; i++)
	  {
	    first = first_element_dof[i];
	    for (j = 0; j < cnt[nv+ned+i]; j++)
	      table[nv+ned+i][j] = first+j;
	  }


	/*
	// VE + I

	int i, j, first;
	ARRAY<int> cnt(nv+nel);
	cnt = 0;
	for (i = 0; i < nv; i++)
	  cnt[i] = 1;
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    int ndof = first_edge_dof[i+1]-first_edge_dof[i];
	    ma.GetEdgePNums (i, v1, v2);
	    cnt[v1] += ndof;
	    cnt[v2] += ndof;
	  }
	for (i = 0; i < nel; i++)
	  cnt[nv+i] = first_element_dof[i+1]-first_element_dof[i];
	
	Table<int> & table = *new Table<int> (cnt);
	cnt = 0;
	for (i = 0; i < nv; i++)
	  {
	    table[i][0] = i;
	    cnt[i] = 1;
	  }
	for (i = 0; i < ned; i++)
	  {
	    int v1, v2;
	    first = first_edge_dof[i];
	    int ndof = first_edge_dof[i+1]-first;
	    ma.GetEdgePNums (i, v1, v2);
	    for (j = 0; j < ndof; j++)
	      {
		table[v1][cnt[v1]++] = first+j;
		table[v2][cnt[v2]++] = first+j;
	      }
	  }
	for (i = 0; i < nel; i++)
	  {
	    first = first_element_dof[i];
	    int ndof = first_element_dof[i+1]-first;
	    for (j = 0; j < ndof; j++)
	      table[nv+i][j] = first+j;
	  }
	*/



	//(*testout) << "H1HO-table = " << table << endl;
	return &table;
      }
    return 0;
  }


  ARRAY<int> * 
  H1HighOrderFESpace :: CreateDirectSolverClusters (int type) const
  {

    int i, j, k;

    int nv = ma.GetNV();
    int nd = GetNDof();
    int ne = ma.GetNE();
    ARRAY<int> & clusters = *new ARRAY<int> (GetNDof());
    clusters = 0;

    /*  // all vertices in global space
    
    for (i = 0; i < nv; i++)
      clusters[i] = 1;

    // all edges in direct solver cluster

    for (i = first_edge_dof[0]; i < first_edge_dof.Last(); i++)
      clusters[i] = 1;
    return &clusters;

    */
   
    // All Vertical Edges in one Cluster for Hex and Prism (-> 2d Problems !) 
    ARRAY<int> ednums, edorient,fnums, forient;

    //ARRAY<int> & clusters = *new ARRAY<int> (nd);
    //clusters = 0;


    
    for (i = 0; i < ne; i++)
      {
	if (ma.GetElType(i) == ET_PRISM)
	  {
	    ma.GetElEdges (i, ednums, edorient);
	    for (j = 6; j < 9; j++)  //vertical Edges 
	      { 
		int first = first_edge_dof[ednums[j]];
		int next = first_edge_dof[ednums[j]+1];
		for (k = first; k < next; k++)
		  clusters[k] = 2;
	      }
	    
	    ma.GetElFaces(i,fnums,forient); // vertical faces 
	      for(j=2;j<5;j++) 
	      {
	
		int first = first_face_dof[fnums[j]]; 
		int next = first_face_dof[fnums[j]+1]; 
		
		for (k=first; k < next; k++) 
		  clusters[k]=0; 
		
		//	int p = order_face[fnums[j]];
		// for(k=2*(p+1)*(p+1);k<next;k++)
		//   clusters[k]=3;  



	      }


	    
	  }

	else if (ma.GetElType(i) == ET_HEX)
	  {
	    ma.GetElEdges (i, ednums, edorient);
	    for (j = 8; j < 12; j++)
	      {
		int first = first_edge_dof[ednums[j]];
		int next = first_edge_dof[ednums[j]+1];
		for (k = first; k < next; k++)
		  clusters[k] = 2;
	      }
	     ma.GetElFaces(i,fnums,forient); // vertical faces 
	      for(j=2;j<6;j++) 
	      {
	
		int first = first_face_dof[fnums[j]]; 
		int next = first_face_dof[fnums[j]+1]; 
		
		for (k=first; k < next; k++) 
		  clusters[k]=0; 
	      }

	  } 
      }


    
    for(i=0; directsolverclustered.Size() > 0 && i<ne; i++)
      {
	if(directsolverclustered[ma.GetElIndex(i)])
	  {
	    GetDofNrs(i,ednums);
	    for(k=0; k<ednums.Size(); k++)
	      {
		clusters[ednums[k]] = 3;
	      }
	  }
      }
       
    // (*testout) << "direct clusters = " << endl << clusters << endl;

    return &clusters;
    
  }




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