/**
   High Order Finite Element Space for H(Curl) 
*/
#include <comp.hpp>
#include <fem.hpp> 

namespace ngcomp 
{
  using namespace ngcomp; 

  HCurlHighOrderFESpace ::  
  HCurlHighOrderFESpace (const MeshAccess & ama, const Flags & aflags)
    : FESpace (ama, aflags),
      flags(aflags)
  {
    gradientdomains.SetSize (ma.GetNDomains());
    gradientdomains.Set();
    
    fn=int(flags.GetNumFlag("face",1));
    
    if (flags.NumListFlagDefined("gradientdomains"))
      {
	const ARRAY<double> & graddomains = flags.GetNumListFlag ("gradientdomains");
	// cout << " graddomains " << graddomains << endl; 
	for (int i = 0; i < gradientdomains.Size(); i++)
	  if (!graddomains[i])
	    gradientdomains.Clear(i);

	cout << endl << " ******* gradientdomains" << gradientdomains << endl;
      }
    
    if(flags.GetDefineFlag("nograds"))
      gradientdomains.Clear(); 
    
    low_order_space = new  NedelecFESpace (ma, 0, dimension, iscomplex);

    rel_order = int (flags.GetNumFlag ("relorder", 0));
    smoother = int (flags.GetNumFlag("smoothing",8));

    

    /*    if (flags.GetDefineFlag ("optext"))
	  {
	  trig = new HCurlHighOrderTrig<TrigExtensionMin> (order);
	  tet =  new HCurlHighOrderTet<TrigExtensionMin> (order);
	  prism =  new HCurlHighOrderPrism<TrigExtensionMin> (order);
	  }
	  else
    */

    trig = new HCurlHighOrderTrig<IntegratedLegendreMonomialExt> (order);
    tet =  new HCurlHighOrderTet<IntegratedLegendreMonomialExt> (order);
    prism =  new HCurlHighOrderPrism<IntegratedLegendreMonomialExt> (order);
    quad = new HCurlHighOrderQuad<IntegratedLegendreMonomialExt> (order);
    hex = new HCurlHighOrderHex<IntegratedLegendreMonomialExt> (order);
    segm = new HCurlHighOrderSegm<IntegratedLegendreMonomialExt> (order); 
    

    // Evaluator for shape tester 
    static ConstantCoefficientFunction one(1);
    if (ma.GetDimension() == 2)
      {
	ARRAY<CoefficientFunction*> coeffs(1);
	coeffs[0] = &one;
	evaluator = GetIntegrators().CreateBFI("massedge", 2, coeffs);
      }
    else if(ma.GetDimension() == 3) 
      {
	ARRAY<CoefficientFunction*> coeffs(1); 
	coeffs[0] = &one;
	evaluator = GetIntegrators().CreateBFI("massedge",3,coeffs); 
	boundary_evaluator = GetIntegrators().CreateBFI("robinedge",3,coeffs); 
	
      }

  }
  
  HCurlHighOrderFESpace :: ~HCurlHighOrderFESpace ()
  {
    ;
  }
  
  FESpace * HCurlHighOrderFESpace :: 
  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 NedelecFESpace (ma,0, dim, iscomplex);      // Space with MG   
    else       
    */
    return new HCurlHighOrderFESpace (ma, flags);
  }
  
  void HCurlHighOrderFESpace :: Update()
  {
    int i, j; 
    
    if (order < 0) 
      throw Exception("HCurlHighOrderFESpace::Update() order < 0 !" ) ;
    
    low_order_space -> Update();
    
    nv = ma.GetNV();
    nel = ma.GetNE();
    
    switch(ma.GetDimension())
      {
      case 2: 
	ned = ma.GetNEdges();
	nfa = 0;  
	break; 
      case 3: 
	ned = ma.GetNEdges();
	nfa = ma.GetNFaces();
	break; 
      }
    
    order_edge.SetSize (ned);   
    order_face.SetSize (nfa); 
    order_inner.SetSize (nel);
    usegrad_edge.SetSize (ned);                                
    usegrad_face.SetSize (nfa); 
    usegrad_cell.SetSize (nel);
     
    order_edge = order;
    order_face = order; 
    order_inner = order;
    // order_inner = 0;
   
    usegrad_edge = 0;                                
    usegrad_face = 0; 
    usegrad_cell = 0; 

    ARRAY<int> eledges, elfaces;
  
     
    for(i=0; i< nel; i++) 
      if(gradientdomains[ma.GetElIndex(i)]) 
	{
	  ma.GetElEdges(i,eledges);
	  	  
	  for(j=0;j<eledges.Size();j++)
	    usegrad_edge[eledges[j]]=1;
	  
	  if(ma.GetDimension()==3)
	    {
	      ma.GetElFaces(i,elfaces);
	      
	      for(j=0;j<elfaces.Size();j++) 
		usegrad_face[elfaces[j]]=1; 
	    }
	  usegrad_cell[i] = 1; 
	}

    /*  if (ma.GetDimension() == 3)
      for (i = 0; i < nel; i++)
	{	  
	  int elorder = ma.GetElOrder(i) - 1 + rel_order; // start with 0
	  //	  cout << "GetElOrder" << ma.GetElOrder(i) << " elorder " << elorder << endl; 
	  elorder = order; 

	  ma.GetElEdges (i, eledges);
	  ma.GetElFaces (i, elfaces);
	  
	  for (j = 0; j < eledges.Size(); j++)
	    if (elorder > order_edge[eledges[j]])
	      order_edge[eledges[j]] = elorder;
	  
	  for (j = 0; j < elfaces.Size(); j++)
	    if (elorder > order_face[elfaces[j]])
	      order_face[elfaces[j]] = elorder;

	  if (elorder > order_inner[i])
	    order_inner[i] = elorder;
	    } */ 


    /*
    // 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]] = 0;
	  order_face[face] = 0;
      }
    */

    /*
    (*testout) << " order = " << order << endl; 
    (*testout) << "order_edge = " << order_edge << endl;
    (*testout) << "order_face = " << order_face << endl;
    (*testout) << "order_inner = " << order_inner << endl;
    */

    ARRAY<int> pnums; 
    
    ndof = ned; // Nedelec (order = 0) !!    
    first_edge_dof.SetSize(ned+1); // last element = ndof;  
    for (i = 0; i < ned; i++)
      {
	first_edge_dof[i] = ndof;
	ndof += usegrad_edge[i]*order_edge[i];
      }
        	
    first_edge_dof[ned] = ndof;
    
    first_face_dof.SetSize(nfa+1);
    int inci =0;
   
    for (i=0; i< nfa; i++) 
      { 
	if(order_face[i]>0) 
	  {
	    ma.GetFacePNums(i,pnums);  
	    int p = order_face[i]; 
	    switch(pnums.Size())
	      {
	      case 3: //Triangle   
		inci = ((usegrad_face[i]+1)*p + 2)*(p-1)/2;    
		break; 
	      case 4: //Quad   
		// inci= 2*(p+1)*p;
		inci = ((usegrad_face[i]+1)*p + 2) * p;  
		break; 
	      }
	  }
	else 
	  inci = 0;
	first_face_dof[i] = ndof; 
	ndof+= inci; 
	
      } 
    first_face_dof[nfa] = ndof;     
    

    first_inner_dof.SetSize(nel+1);
    for (i=0; i< nel; i++)
      {
	int p = order_inner[i];
	switch(ma.GetElType(i)) 
	  {
	  case ET_TRIG:
	    inci= ((usegrad_cell[i]+1)*p + 2) * (p-1) /2;;
	    break; 
	  case ET_QUAD: 
	    inci= ((usegrad_cell[i]+1) * p + 2) * p; 
	    break; 
	  case ET_TET: 
	    inci = ((usegrad_cell[i] + 2) *  p + 3) * (p-2) * (p-1) / 6; 
	    break; 
	  case ET_PRISM:
	    inci = ((usegrad_cell[i]+2)* p + 3) * p * (p-1) / 2; 
	    break; 
	  case ET_HEX:
	    inci =  ((usegrad_cell[i] + 2 )* p + 3) * p * p;
	    break; 
	  default:  
	    inci = 0;  
	    break; 
	  }

	if (inci < 0 || p <= 0) inci = 0;
	first_inner_dof[i] = ndof; 
	ndof+= inci;
      }
    first_inner_dof[nel] = ndof;     
    
    /*
      (*testout) << " HCURL " << endl; 
      (*testout) << "usegrad_edge " << usegrad_edge << endl; 
      (*testout) << "usegrad_face " << usegrad_face << endl; 
      (*testout) << "usegrad_cell " << usegrad_cell << endl; 
      (*testout) << "ndof  = " << ndof << endl; 
      (*testout) << "first_edge   = " << first_edge_dof << endl;
      (*testout) << "first_face   = " << first_face_dof << endl;
      (*testout) << "first_inner  = " << first_inner_dof << endl;
      (*testout) << " order = " << order << endl; 
      (*testout) << "order_edge HC = " << order_edge << endl;
      (*testout) << "order_face HC = " << order_face << endl;
      (*testout) << "order_inner HC = " << order_inner << endl;
    */
  }
  
  const FiniteElement & HCurlHighOrderFESpace :: 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_TRIG:
	fe = trig; break;
      case ET_QUAD:
	fe = quad; break;
      case ET_HEX: 
	fe = hex; break; 
      default:
	fe = 0;
      }
    
    if (!fe)
      {
	stringstream str;
	str << "HCurlHighOrderFESpace " << GetClassName() 
	    << ", undefined eltype " 
	    << ElementTopology::GetElementName(ma.GetElType(elnr))
	    << ", order = " << order << endl;
	throw Exception (str.str());
      }

    ArrayMem<int,12> vnums; // calls GetElPNums -> max 12 for PRISM12
    ma.GetElVertices(elnr, vnums);

    if(ma.GetDimension() == 2) 
      {	
	HCurlHighOrderFiniteElement<2> * hofe = 
	  dynamic_cast<HCurlHighOrderFiniteElement<2>*> (fe);
	
    	ArrayMem<int, 4> ednums, ord_edge, ug_edge;
	
	ma.GetElEdges(elnr, ednums);
	
	ord_edge.SetSize (ednums.Size());
	ug_edge.SetSize (ednums.Size()); 
	
	for (int j = 0; j < ednums.Size(); j++)
	  {
	    ord_edge[j] = order_edge[ednums[j]];
	    ug_edge[j]  = usegrad_edge[ednums[j]]; 
	  }

	hofe -> SetVertexNumbers (vnums);
	hofe -> SetOrderEdge (ord_edge);
	hofe -> SetOrderInner (order_inner[elnr]);

	hofe -> SetUsegradEdge (ug_edge); 
	hofe -> SetUsegradCell (usegrad_cell[elnr]); 

	hofe -> ComputeNDof();
      }   
    else if (ma.GetDimension() == 3) 
      {
	HCurlHighOrderFiniteElement<3> * hofe = 
	  dynamic_cast<HCurlHighOrderFiniteElement<3>*> (fe);

	ArrayMem<int, 12> ednums, ord_edge, ug_edge;
	ArrayMem<int, 6> fanums, ord_face, ug_face;
	
	ma.GetElEdges(elnr, ednums);
	ma.GetElFaces(elnr, fanums);
	
	ord_edge.SetSize (ednums.Size());
	ord_face.SetSize (fanums.Size());
	ug_edge.SetSize (ednums.Size()); 
	ug_face.SetSize (fanums.Size()); 
	
	for (int j = 0; j < ednums.Size(); j++)
	  {
	    ord_edge[j] = order_edge[ednums[j]];
	    ug_edge[j] = usegrad_edge[ednums[j]]; 
	  }
	for (int j = 0; j < fanums.Size(); j++)
	  {
	    ord_face[j] = order_face[fanums[j]];
	    ug_face[j] = usegrad_face[fanums[j]];
	  }

	hofe -> SetVertexNumbers (vnums);
	hofe -> SetOrderEdge (ord_edge);
	hofe -> SetOrderFace (ord_face);
	hofe -> SetOrderInner (order_inner[elnr]);

	hofe -> SetUsegradEdge (ug_edge); 
	hofe -> SetUsegradFace (ug_face); 
	hofe -> SetUsegradCell (usegrad_cell[elnr]); 

	hofe -> ComputeNDof();
      }

    return *fe;
  }
 
  const FiniteElement & HCurlHighOrderFESpace :: GetSFE (int selnr, LocalHeap & lh) const
  {
    int i, j;

    FiniteElement * fe = 0;
    
    switch (ma.GetSElType(selnr))
      {
      case ET_SEGM: 
	fe = segm; break; 
      case ET_TRIG:
	fe = trig; break;
      case ET_QUAD:
	fe = quad; break;
      default:
	fe = 0;
      }
    
    if (!fe)
      {
	stringstream str;
	str << "HCurlHighOrderFESpace " << GetClassName() 
	    << ", undefined eltype " 
	    << ElementTopology::GetElementName(ma.GetSElType(selnr))
	    << ", order = " << order << endl;
	throw Exception (str.str());
      }

    ArrayMem<int, 8> vnums;
    ArrayMem<int, 4> ednums, ord_edge, ug_edge; 
    int ord_face, ug_face; 

    ma.GetSElVertices(selnr, vnums);

    if(fe == segm) 
      {
	HCurlHighOrderFiniteElement<1> * hofe =
	  dynamic_cast<HCurlHighOrderFiniteElement<1>*> (fe);
	
	hofe -> SetVertexNumbers (vnums);
	
	ma.GetSElEdges(selnr, ednums);
	hofe -> SetOrderInner (order_edge[ednums[0]]);
	hofe -> SetUsegradCell (usegrad_edge[ednums[0]]);
	hofe -> ComputeNDof();
      }
    else 
      {     
	HCurlHighOrderFiniteElement<2> * hofe =
	  dynamic_cast<HCurlHighOrderFiniteElement<2>*> (fe);
      
	ma.GetSElEdges(selnr, ednums);
	ord_face = order_face[ma.GetSElFace(selnr)];
	ug_face = usegrad_face[ma.GetSElFace(selnr)];
	
	ord_edge.SetSize (ednums.Size());
	ug_edge.SetSize (ednums.Size()); 
	
	for (j = 0; j < ednums.Size(); j++)
	  {
	    ord_edge[j] = order_edge[ednums[j]]; 
	    ug_edge[j] = usegrad_edge[ednums[j]];
	  }
	
	hofe -> SetVertexNumbers (vnums);
	hofe -> SetOrderEdge (ord_edge);
	hofe -> SetOrderInner (ord_face);
      
	hofe -> SetUsegradEdge(ug_edge); 
	hofe -> SetUsegradCell(ug_face); 
      
	hofe -> ComputeNDof();
	
      }
    
    return *fe;
  }
 


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

  void HCurlHighOrderFESpace :: GetDofNrs (int elnr, ARRAY<int> & dnums) const
  {
    // Save phase-functions the way 
    // (1*e1),.. (1*e_ne)  
    // (p_t)*tang_e1, ... p_t*tang_ne
    // p_n * el1, p_i * el1, ... , p_n * nel_n , p_i *el_nel 

    ARRAY<int> vnums, ednums, fnums;
    int i, j;
    int first,next; 

    ma.GetElEdges (elnr, ednums);
    if(ma.GetDimension() == 3)
      ma.GetElFaces (elnr, fnums); 
    else 
      fnums.SetSize(0); 
    dnums.SetSize(0);

 
    if(order < 0) 
      throw Exception(" HCurlHighOrderFESpace :: GetDofNrs() order < 0 "); 

    //Nedelec0
    for (i = 0; i < ednums.Size(); i++) 
      dnums.Append (ednums[i]);
    
    //edges_tang
    for (i = 0; i < ednums.Size(); i++) 
      {
	first = first_edge_dof[ednums[i]];
	next = first_edge_dof[ednums[i]+1]; 
	
	for (j = first; j <next; j++)
	  dnums.Append (j);
      }
    
    
    // faces 
    for(i=0; i<fnums.Size(); i++)     
      { 
	first = first_face_dof[fnums[i]]; 
	next = first_face_dof[fnums[i]+1];
	for(j=first ; j<next; j++) 
	  dnums.Append(j); 
      } 
    
    //inner 
    first = first_inner_dof[elnr];
    next = first_inner_dof[elnr+1]; 

    for(j=first; j<next; j++) 
      dnums.Append(j);
  }
 

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

    // Save phase-functions the way 
    // (1*e1),.. (1*e_ne)  
    // (p_t)*tang_e1, ... p_t*tang_ne
    // p_n * el1, p_i * el1, ... , p_n * nel_n , p_i *el_nel 

    ARRAY<int> vnums, ednums, fnums;
    int i, j;
    int first,next; 

    ma.GetElEdges (elnr, ednums);
    if(ma.GetDimension() == 3)
      ma.GetElFaces (elnr, fnums); 
    else 
      fnums.SetSize(0); 
    dnums.SetSize(0);

 
    if(order < 0) 
      throw Exception(" HCurlHighOrderFESpace :: GetDofNrs() order < 0 "); 

    //Nedelec0
    for (i = 0; i < ednums.Size(); i++) 
      dnums.Append (ednums[i]);
    
    //edges_tang
    for (i = 0; i < ednums.Size(); i++) 
      {
	first = first_edge_dof[ednums[i]];
	next = first_edge_dof[ednums[i]+1]; 
	
	for (j = first; j <next; j++)
	  dnums.Append (j);
      }
    
    
    // faces 
    if (ma.GetDimension() == 3)
      for(i=0; i<fnums.Size(); i++)     
	{ 
	  first = first_face_dof[fnums[i]]; 
	  next = first_face_dof[fnums[i]+1];
	  for(j=first ; j<next; j++) 
	    dnums.Append(j); 
	} 
  }
 


  void HCurlHighOrderFESpace :: GetSDofNrs (int selnr, ARRAY<int> & dnums) const
  {
    // Save shape-functions the way 
    // (1*e1),.. (1*e_ne)  
    // (p_t)*tang_e1, ... p_t*tang_ne
    // p_n * el1, p_i * el1, ... , p_n * nel_n , p_i *el_nel 

    ARRAY<int> vnums, ednums;
    int fnum; 
    int i, j;
    int first,next; 

    ma.GetSElEdges (selnr, ednums);
    fnum = ma.GetSElFace (selnr); 

    dnums.SetSize(0);

    if(order <0) throw (" HCurlHighOrderFESpace :: GetSDofNrs() order < 0 "); 
    //Nedelec
    for (i = 0; i < ednums.Size(); i++) 
      dnums.Append (ednums[i]);

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

    if(first_face_dof.Size()>1)		       
      {
	// inner = normal + bubbles 
	first = first_face_dof[fnum];
	next = first_face_dof[fnum+1]; 
	for(j=first; j<next; j++) 
	  dnums.Append(j);
      }
    
    // (*testout) << "surf-dnums = " << endl << dnums << endl;
  }
 
  
  Table<int> * HCurlHighOrderFESpace :: 
  CreateSmoothingBlocks (int type) const
  {
    int i, j, k, first;
    int ncnt; 
    
    int ii; 

    //    int SmoothingType=7; 
    int SmoothingType = smoother; 

    // 1: ned-block, ho-edes, inner 
    // 2: neds, ho-edges, inner 
    // 3: ned+ho-edges, inner
    // 4: ho-edges, inner (no ned Block)
    // 5: vertices with lo-edges, ho-edges, inner
    // 6: ho-edges, faces+inner 
    // 7: Clustering: ho-edges, faces, inner
    // 8: vertices with lo+ho edges, faces, inner 
    // 9: vertex augmented
   
    cout << "SmoothingType " << SmoothingType << endl; 
    
    ARRAY<int> vnums; 
    ARRAY<int> orient; 
    ARRAY<int> ednums, fanums, faorient;

    for(i=0; i< nel ; i++)                                                                                                                                                      ma.GetElVertices (i, vnums);
      
    for(i=0; i< nel; i++)  
      ma.GetElEdges (i,vnums,orient);
	      
    int pn1,pn2; 
    for(i=0; i< ned; i++) 
      ma.GetEdgePNums (i,pn1,pn2);

    switch(SmoothingType) 
      {
      case 1: 
	// 1 RT Block 
	ncnt = 1 + first_edge_dof.Size() -1 + first_inner_dof.Size() -1; 	    
	break; 
	
      case 2: 
	// ned RT Blocks 
	ncnt = ned + first_edge_dof.Size() -1 + first_inner_dof.Size() -1; 	    
	break; 
	
      case 3: 
	// ned RT+ edges_normal Blocks   
	ncnt = ned + first_inner_dof.Size() -1; 	  
	break; 
      case 4: 
	// no  RT Block 
	ncnt = first_edge_dof.Size() -1 + first_inner_dof.Size() -1; 	    
	break; 
	
      case 5: 
	// RT in divergence blocks 
	ncnt = nv + first_edge_dof.Size() -1 + first_inner_dof.Size() -1; 
	break; 
      case 6:
	ncnt = ned + nfa;
	break;
      case 7: 
	ncnt = nv + ned + nfa + nel;
	break;
      case 8: 
	ncnt = nv + ned + nfa + nel;
	break;
      case 9: 
	ncnt = nv + ned + nfa + nel;
	break;
      case 10:
	ncnt = nv + ned + nfa + nel; 
      default: 
	cout << " Error in HCurlHOFESpace :: CreateSmoothingBlocks no SmoothingType " << endl; 
	ncnt =0; 
      } 

    ARRAY<int> cnt(ncnt); 
    cnt = 0;
    
    cout << " ncnt " << ncnt << endl; 

    switch (SmoothingType)
      {
      case 1: 
	ii = 0;  
	cnt[ii] = ned; //RT
	ii++;
	for(i=0; i< first_edge_dof.Size()-1;i++ )
	  cnt[ii++] = first_edge_dof[i+1] - first_edge_dof[i]; 
	for(i=0; i< first_inner_dof.Size()-1;i++ )
	  cnt[ii++] = first_inner_dof[i+1] - first_inner_dof[i]; 	
	break; 
      case 2: 
	ii = 0;  
	for(i=0; i< ned; i++)  cnt[ii++] = 1; 
	for(i=0; i< first_edge_dof.Size()-1;i++ )
	  cnt[ii++] = first_edge_dof[i+1] - first_edge_dof[i]; 
	for(i=0; i< first_inner_dof.Size()-1;i++ )
	  cnt[ii++] = first_inner_dof[i+1] - first_inner_dof[i]; 
	break; 
      case 3: 
	ii = 0;  
	for(i=0; i< ned; i++)  cnt[i] = 1; 
	for(i=0; i< first_edge_dof.Size()-1;i++ )
	  cnt[i] += first_edge_dof[i+1] - first_edge_dof[i];
	ii=ned;
	for(i=0; i< first_inner_dof.Size()-1;i++)
	  cnt[ii++] = first_inner_dof[i+1] - first_inner_dof[i]; 
	break; 
      case 4: 
	ii = 0;  
	for(i=0; i< first_edge_dof.Size()-1;i++ )
	  cnt[ii++] = first_edge_dof[i+1] - first_edge_dof[i]; 
	for(i=0; i< first_inner_dof.Size()-1;i++ )
	  cnt[ii++] = first_inner_dof[i+1] - first_inner_dof[i]; 	
	break; 	
      case 5: 
	int v1, v2; 
	for(i=0; i< nv; i++) cnt[i]=0;  
	for(i=0; i< ned; i++) 
	  {
	    ma.GetEdgePNums (i, v1, v2);
	  
	    cnt[v1]++; 
	    cnt[v2]++;
	   
	  }
	ii = nv; 
	for(i=0; i< first_edge_dof.Size()-1;i++ )
	  cnt[ii++] = first_edge_dof[i+1] - first_edge_dof[i]; 
	for(i=0; i< first_inner_dof.Size()-1;i++ )
	  cnt[ii++] = first_inner_dof[i+1] - first_inner_dof[i]; 	
	break; 	
      case 6:
	{
	  for (i = 0; i < ned; i++)
	    cnt[i] = first_edge_dof[i+1] - first_edge_dof[i];
	  for (i = 0; i < nfa; i++)
	    cnt[ned+i] = first_face_dof[i+1] - first_face_dof[i];
	  for (i = 0; i < nel; i++)
	    {	   
	      if(ma.GetDimension() == 3) ma.GetElFaces (i, fanums, faorient);
	      else fanums.SetSize(0); 
	      for (j = 0; j < fanums.Size(); j++)
		cnt[ned+fanums[j]] += first_inner_dof[i+1]-first_inner_dof[i];
	    }
	  break;
	}
      case 7:
	{
	  /*  
	  for (i = 0; i < ned; i++)
	    cnt[i] = first_edge_dof[i+1] - first_edge_dof[i];
	  for (i = 0; i < nfa; i++)
	    cnt[ned+i] = first_face_dof[i+1] - first_face_dof[i];
	  for (i = 0; i < nel; i++)
	    cnt[ned+nfa+i] = first_inner_dof[i+1] - first_inner_dof[i];
	  break;
	  */
	  for (i = 0; i < ned; i++)
	    cnt[ma.GetClusterRepEdge(i)] +=
	      first_edge_dof[i+1] - first_edge_dof[i];
	  for (i = 0; i < nfa; i++)
	    cnt[ma.GetClusterRepFace(i)] += 
	      first_face_dof[i+1] - first_face_dof[i];
	  if (!eliminate_internal)
	    for (i = 0; i < nel; i++)
	      cnt[ma.GetClusterRepElement(i)] += 
		first_inner_dof[i+1] - first_inner_dof[i];
	  break;
	  

	}
      case 8:
	{
	  for (i = 0; i < ned; i++)
	    {
	      ma.GetEdgePNums (i,pn1,pn2);
	      cnt[pn1] += 1 + first_edge_dof[i+1] - first_edge_dof[i];
	      cnt[pn2] += 1 + 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_inner_dof[i+1] - first_inner_dof[i];
	  else
	    for (i = 0; i < nel; i++)
	      cnt[nv+ned+nfa+i] = 0;
	  break;
	}
      case 9:
	{
	  for (i = 0; i < nv; i++)
	    cnt[i] = 1;
	  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_inner_dof[i+1] - first_inner_dof[i];

	      ma.GetElFaces (i,fanums,orient);
	      for (k = 0; k < fanums.Size(); k++)
		cnt[nv+ned+fanums[k]] += first_inner_dof[i+1] - first_inner_dof[i];
	    }
	  break;
	}
      }
    
 
    Table<int> & table = *new Table<int> (cnt); 
    
    switch(SmoothingType) 
      {
      case 1: 
	ii=0; 
	for(i=0; i<ned; i++) table[ii][i] = i; //RT 
	ii++;
	for(i=0; i<first_edge_dof.Size()-1; i++) 
	  {
	    first = first_edge_dof[i]; 
	    for(j=0; j< cnt[ii]; j++)
	      table[ii][j] = j+first; 
	    ii++; 
	  } 
	for(i=0; i<first_inner_dof.Size()-1; i++) 
	  { 
	    first = first_inner_dof[i]; 
	    for(j=0; j< cnt[ii]; j++) 
	      table[ii][j] = j+ first; 
	    ii++;
	  } 	
	break; 
	
      case 2: 
	ii=0; 
	for(i=0; i<ned; i++) table[ii++][0] = i; //RT 
	for(i=0; i<first_edge_dof.Size() -1;i++) 
	  {
	    first = first_edge_dof[i]; 
	    for(j=0; j< cnt[ii]; j++)
	      table[ii][j] = j+first; 
	    ii++; 
	  } 
	for(i=0; i<first_inner_dof.Size()-1;i++) 
	  { 
	    first = first_inner_dof[i]; 
	    for(j=0; j< cnt[ii]; j++) 
	      table[ii][j] = j+ first; 
	    ii++;
	  } 
	break; 

      case 3: 
	ii=0; 
	for(i=0; i<ned; i++) 
	  { 
	    table[ii++][0] = i; //RT 
	  }
	for(i=0; i<first_edge_dof.Size()-1;i++) 
	  {
	    first = first_edge_dof[i]; 
	    for(j=1; j< cnt[i]; j++)
	      table[i][j] = j+first-1;  
	  } 
	for(i=0; i<first_inner_dof.Size()-1;i++) 
	  { 
	    first = first_inner_dof[i]; 
	    for(j=0; j< cnt[ii]; j++) 
	      table[ii][j] = j+ first; 
	    ii++;
	  } 
	break; 
      case 4: 
	ii=0; 
	for(i=0; i<first_edge_dof.Size()-1; i++) 
	  {
	    first = first_edge_dof[i]; 
	    for(j=0; j< cnt[ii]; j++)
	      table[ii][j] = j+first; 
	    ii++; 
	  } 
	for(i=0; i<first_inner_dof.Size()-1; i++) 
	  { 
	    first = first_inner_dof[i]; 
	    for(j=0; j< cnt[ii]; j++) 
	      table[ii][j] = j+ first; 
	    ii++;
	  } 	
	break;
      case 5:
	int v1,v2; 
	for(i=0; i< ned; i++) 
	  {
	    ma.GetEdgePNums (i, v1, v2);
	    table[v1][--cnt[v1]] = i; 
	    table[v2][--cnt[v2]] = i; 
	  }
	ii = nv; 
	for(i=0; i<first_edge_dof.Size()-1; i++) 
	  {
	    first = first_edge_dof[i]; 
	    for(j=0; j< cnt[ii]; j++)
	      table[ii][j] = j+first; 
	    ii++; 
	  } 
	for(i=0; i<first_inner_dof.Size()-1; i++) 
	  { 
	    first = first_inner_dof[i]; 
	    for(j=0; j< cnt[ii]; j++) 
	      table[ii][j] = j+ first; 
	    ii++;
	  } 
	break; 
      case 6:
	{
	  cnt = 0; 
	  for (i = 0; i < ned; i++)
	    for (j = first_edge_dof[i]; j < first_edge_dof[i+1]; j++)
	      table[i][cnt[i]++] = j;
	  for (i = 0; i < nfa; i++)
	    for (j = first_face_dof[i]; j < first_face_dof[i+1]; j++)
	      table[ned+i][cnt[ned+i]++] = j;
	  for (i = 0; i < nel; i++)
	    {
	      if(ma.GetDimension() == 3)  ma.GetElFaces (i, fanums, faorient);
	      else fanums.SetSize(0); 
	      for (j = 0; j < fanums.Size(); j++)
		for (k = first_inner_dof[i]; k < first_inner_dof[i+1]; k++)
		  table[ned+fanums[j]][cnt[ned+fanums[j]]++] = k;
	    }
	  break;
	}
      case 7:
	{
	  cnt = 0; 
	  for (i = 0; i < ned; i++)
	    {
	      int nr = ma.GetClusterRepEdge (i);
	      for (j = first_edge_dof[i]; j < first_edge_dof[i+1]; j++)
		table[nr][cnt[nr]++] = j;
	    }
	  for (i = 0; i < nfa; i++)
	    {
	      int nr = ma.GetClusterRepFace (i);
	      for (j = first_face_dof[i]; j < first_face_dof[i+1]; j++)
		table[nr][cnt[nr]++] = j;
	    }
	  for (i = 0; i < nel; i++)
	    {
	      int nr = ma.GetClusterRepElement (i);
	      for (j = first_inner_dof[i]; j < first_inner_dof[i+1]; j++)
		table[nr][cnt[nr]++] = j;
	    }

	  /*	  if (!eliminate_internal)
	    for (i = 0; i < nel; i++)
	      {
		int nr = ma.GetClusterRepElement (i);
		for (j = first_inner_dof[i]; j < first_inner_dof[i+1]; j++)
		  table[nr][cnt[nr]++] = j;
	      }

	  */
	  
	  /*
	  cnt = 0;
	  for (i = 0; i < ned; i++)
	    for (j = first_edge_dof[i]; j < first_edge_dof[i+1]; j++)
	      table[i][cnt[i]++] = j;
	  for (i = 0; i < nfa; i++)
	    for (j = first_face_dof[i]; j < first_face_dof[i+1]; j++)
	      table[ned+i][cnt[ned+i]++] = j;
	  for (i = 0; i < nel; i++)
	    for (j = first_inner_dof[i]; j < first_inner_dof[i+1]; j++)
	      table[ned+nfa+i][cnt[ned+nfa+i]++] = j;
	  */ 
	  break;
	}
      case 8:
	{
	  cnt = 0;
	  for (i = 0; i < ned; i++)
	    {
	      ma.GetEdgePNums (i,pn1,pn2);	      
	      table[pn1][cnt[pn1]++] = i;
	      table[pn2][cnt[pn2]++] = i;
	      for (j = first_edge_dof[i]; j < first_edge_dof[i+1]; j++)
		{
		  table[pn1][cnt[pn1]++] = j;
		  table[pn2][cnt[pn2]++] = j;
		}
	    }
	  for (i = 0; i < nfa; i++)
	    for (j = first_face_dof[i]; j < first_face_dof[i+1]; j++)
	      table[nv+ned+i][cnt[nv+ned+i]++] = j;
	  if (!eliminate_internal)
	    for (i = 0; i < nel; i++)
	      for (j = first_inner_dof[i]; j < first_inner_dof[i+1]; j++)
		table[nv+ned+nfa+i][cnt[nv+ned+nfa+i]++] = j;
	  break;
	}
      case 9:
	{
	  cnt = 0;
	  for (i = 0; i < nv; i++)
	    table[i][cnt[i]++] = nv+i;
	  for (i = 0; i < ned; i++)
	    for (j = first_edge_dof[i]; j < first_edge_dof[i+1]; j++)
	      table[nv+i][cnt[nv+i]++] = j;
	  for (i = 0; i < nfa; i++)
	    for (j = first_face_dof[i]; 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_inner_dof[i]; j < first_inner_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_inner_dof[i]; j < first_inner_dof[i+1]; j++)
		    table[dof][cnt[dof]++] = j;
		}
	    }
	 // (*testout) << "table = " << table << endl;
	  break;
	}
      }
  //(*testout) << "table = " << table << endl;	
     
    return & table;
  }
  
  ARRAY<int> * 
  HCurlHighOrderFESpace :: 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 Vertical Edges in one Cluster for Hex and Prism + const_x*poly faces + const_y*polx faces (-> 2d Problems !) 
    
    ARRAY<int> ednums, edorient,fnums,forient,pnums;

   

    //lo 
    for(i=0;i<ma.GetNEdges();i++)
      clusters[i]=0; 

    for (i = 0; i < ne; i++)
      {
	ma.GetElPNums(i,pnums); 
	if (ma.GetElType(i) == ET_PRISM)
	  {
	    ma.GetElEdges (i, ednums, edorient);
	     for (j = 0; j < 6; j++)  //horizontal Edges 
	      { 
		int first = first_edge_dof[ednums[j]];
		int next = first_edge_dof[ednums[j]+1];
		for (k = first; k < next; k++)
		  clusters[k] = 0;
		clusters[ednums[j]]=0; 
	      }


	    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] = 3;
		clusters[ednums[j]]=0; 
	      }
	    ma.GetElFaces (i, fnums, forient); 
	    
	    for (j=2; j<5 ; j++)  // vertical faces
	      {
		
		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;  
		
	
	      }
	     for (j=0; j<2 ; j++)  // horizontal faces
	      {
		
		int first = first_face_dof[fnums[j]]; 
		int next = first_face_dof[fnums[j]+1]; 
		
		for (k=first; k < next; k++) 
		  clusters[k]=0; 
	      }
	  }

	else if (ma.GetElType(i) == ET_HEX)
	  {
	    ma.GetElEdges (i, ednums, edorient);
	   
	    for(j=0;j<8;j++) // horizontal edges
	       {	
		int first = first_edge_dof[ednums[j]];
		int next = first_edge_dof[ednums[j]+1];
		for (k = first; k < next; k++)
		  clusters[k] = 0;
		
		clusters[ednums[j]]=0; 

	      }
	    

	    for (j = 8; j < 12; 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] = 3;
		clusters[ednums[j]]=0; 
	      }

	    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; 
		
		int p = order_face[fnums[j]];
		for(k=2*(p+1)*(p+1);k<next;k++)
		  clusters[k]=3;  



	      }
	      for(j=0;j<2;j++)  //horizontal faces 
	      {
	
		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(k=first_inner_dof[i];k<first_inner_dof[i+1];k++) 
          clusters[k]=0; 
	
      }
    
    //  (*testout) << "direct clusters = " << endl << clusters << endl;
    
    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]] = 4;
	      }
	  }
      }
    
    return &clusters;
    
  }
  


    
  
  /*  BitArray * HCurlHighOrderFESpace :: 
  CreateIntermediatePlanes (int type) const
  {
    int i;
    ARRAY<int> vnums;
 
    //bool has_cluster = 0;
     for (i = 0; i < ned; i++)
      {
	int pi1, pi2;
	ma.GetEdgePNums (i, pi1, pi2);
	
	if (ma.GetClusterRepVertex (pi1) ==
	    ma.GetClusterRepVertex (pi2))
	  has_cluster = 1;
      }

    if (!has_cluster) return 0;
    
   
    BitArray & ba = *new BitArray (GetNDof());
    ba.Clear();

    for (i = 0; i < ned; i++)
      {
	int pi1, pi2;
	ma.GetEdgePNums (i, pi1, pi2);
	
	if (ma.GetClusterRepVertex (pi1) ==
	    ma.GetClusterRepVertex (pi2))
	  {
	    ba.Set (i);
	    for (int l = first_edge_dof[i];
		 l < first_edge_dof[i+1]; l++)
	      ba.Set (l);
	  }
      }
    
    
    for (i = 0; i < nfa; i++)
      {
	ma.GetFacePNums (i, vnums);
	if (vnums.Size() == 4)
	  {
	    for (int l = first_face_dof[i];
		 l < first_face_dof[i+1]; l++)
	      ba.Set (l);
	  }	    
	  }
    return &ba;
    }*/


  SparseMatrix<double> * 
  HCurlHighOrderFESpace :: CreateGradient() const
  {
    // bool bound = bli.BoundaryForm(); 
    Flags flags2 (flags);
    flags2.SetFlag("order", order+1);
    (*testout) << " ***  Create H1HighOrder for Grad " << endl; 
    H1HighOrderFESpace  fesh1(ma, flags2); 
    fesh1.Update();
    (*testout) << " *** H1HighOrderFe Updated " << endl; 

    int nfa;
    
    if(ma.GetDimension()==3) nfa= ma.GetNFaces(); 
    else nfa = 0; 
    
    int ned = ma.GetNEdges(); 
    int nv  = ma.GetNV(); 
    int nel = ma.GetNE(); 
    
    int dim     = fesh1.GetDimension();
    int dimcurl = GetDimension();

    ARRAY<int> dnums_h1l; 
    ARRAY<int> dnums_hcl;
    LocalHeap lh (1000000);
    
    // Matrix Graph for gradient matrix , start (end) position for assembling 
    ARRAY<int> cnts(ndof);

    // for(int i=0; i<ndof; i++) cnts[i] = 0; 
         
    /*   //each Nedelec shape by two vertices 
    for(int i=0; i< ned; i++) 	cnts[i] = 2;

    //hcurl-tangential-edge-shapes by contruction grad(h1_edges) 
    for(int k=first_edge_dof[0]; k < first_edge_dof[ned]; k++) cnts[k] = 1;
    */ 
 
    // vertices - nedelec
    for (int i = 0; i < ned; i++) 
      cnts[i] = 2; 
    
    for(int i=ned; i<ndof;i++) cnts[i]=0;
    
    for(int i=0; i<ned; i++)
      {
	int l = first_edge_dof[i]; 
	for( int k = fesh1.GetFirstEdgeDof(i); k< fesh1.GetFirstEdgeDof(i+1); 
	     k++, l++)
	  cnts[l] = 1; 
      }
    
    for(int i=0; i< nfa; i++) 
      {
	int l= first_face_dof[i]; 
	for (int k = fesh1.GetFirstFaceDof(i); k<fesh1.GetFirstFaceDof(i+1); 
	     k++, l++) 
	  cnts[l] = 1; 
      }
    
    for(int i=0; i<nel; i++)
      {
	int l= first_inner_dof[i]; 
	for (int k = fesh1.GetFirstElementDof(i); k<fesh1.GetFirstElementDof(i+1); 
	     k++, l++) 
	  cnts[l] = 1; 
      }
  
    //sparse matrix with above matrix graph cnts 
    SparseMatrix<double> & grad = *new SparseMatrix<double>(cnts);  
    
    //cout << " grad cnts " << endl; 
    // (*testout) << " grad cnts " << endl; 

    // vertices - nedelec
    for (int i = 0; i < ned; i++) 
      {
	int p1,p2; 
	ma.GetEdgePNums(i,p1,p2); 
	grad.CreatePosition(i,p1); 
	grad.CreatePosition(i,p2); 
	//if(p1 > p2)  // A&C 
	if (p1 < p2) // new 
	  {
	    grad(i,p1) = -1.;
	    grad(i,p2) = 1.; 
	  }
	else
	  {
	    grad(i,p1) = 1.; 
	    grad(i,p2) = -1.; 
	  }
      } 

    for(int i=0; i<ned; i++)
      {
	int l = first_edge_dof[i]; 
	for( int k = fesh1.GetFirstEdgeDof(i); k< fesh1.GetFirstEdgeDof(i+1); 
	     k++, l++)
	  {
	    grad.CreatePosition(l,k); 
	    grad(l,k)=1.; 
	  }
      }
    
    for(int i=0; i< nfa; i++) 
      {
	int l= first_face_dof[i]; 
	for (int k = fesh1.GetFirstFaceDof(i); k<fesh1.GetFirstFaceDof(i+1); 
	     k++, l++) 
	  {
	    grad.CreatePosition(l,k); 
	    grad(l,k)=1.;
	  }
      }
    
     
    for(int i=0; i<nel; i++)
      {
	int l= first_inner_dof[i]; 
	for (int k = fesh1.GetFirstElementDof(i); k<fesh1.GetFirstElementDof(i+1); 
	     k++, l++) 
	  {
	    grad.CreatePosition(l,k); 
	    grad(l,k)=1.;
	  }	
      }
    //(*testout) << " Global grad " << grad << endl; 
    return &grad;
  }



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