/*********************************************************************/
/* File:   meshaccess.cpp                                            */
/* Author: Joachim Schoeberl                                         */
/* Date:   25. Mar. 2000                                             */
/*********************************************************************/

/* 
   Access to fe mesh
*/


#include <comp.hpp>

namespace ngcomp

{
  using namespace ngcomp;


  MeshAccess :: MeshAccess ()
  {
    ;
  }

  MeshAccess :: ~MeshAccess ()
  {
    ;
  }


  int MeshAccess :: GetNDomains () const
  {
    int maxind = -1;
    int ne = GetNE(); 
    for (int i = 0; i < ne; i++)
      { 
	int ind = GetElIndex(i);
	if (ind > maxind) maxind = ind;
      }
    return maxind+1;
  }


  int MeshAccess :: GetNBoundaries () const
  {
    int maxind = -1;
    int nse = GetNSE(); 
    for (int i = 0; i < nse; i++)
      { 
	int ind = GetSElIndex(i);
	if (ind > maxind) maxind = ind;
      }
    return maxind+1;
  }


  
ELEMENT_TYPE MeshAccess :: GetElType (int elnr) const
{
  elnr++;
  switch (Ng_GetElementType (elnr))
    {
    case NG_TRIG:
    case NG_TRIG6:
      {
	return ET_TRIG;
	break;
      }

    case NG_QUAD:
    case NG_QUAD6:
      {
	return ET_QUAD;
	break;
      }

    case NG_TET:
    case NG_TET10:
      {
	return ET_TET;
	break;
      }
    case NG_PRISM:
    case NG_PRISM12:
      {
	return ET_PRISM;
	break;
      }
    case NG_PYRAMID:
      {
	return ET_PYRAMID;
	break;
      }
    case NG_HEX:
      {
	return ET_HEX;
	break;
      }

    default:
      {
	cerr << "MeshAccess::GetElType undefined type "
	     << int(Ng_GetElementType (elnr)) << endl;
      }
    }
  return ET_TET;
}


ELEMENT_TYPE MeshAccess :: GetSElType (int elnr) const
{
  int pi[NG_SURFACE_ELEMENT_MAXPOINTS];
  switch (Ng_GetSurfaceElement (elnr+1, pi))
    {
    case NG_SEGM:
    case NG_SEGM3:
      {
	return ET_SEGM;
	break;
      }

    case NG_TRIG:
    case NG_TRIG6:
      {
	return ET_TRIG;
	break;
      }

    case NG_QUAD:
    case NG_QUAD6:
      {
	return ET_QUAD;
	break;
      }
    default:
      {
	throw Exception ("GetSElType: Unhandled element type");
      }
    }


  stringstream str;
  str << "MeshAccess::GetSElType undefined type "
      << int(Ng_GetSurfaceElement (elnr, pi)) << endl;
  throw Exception (str.str());
}

  
void MeshAccess :: GetElPNums (int elnr, ARRAY<int> & pnums) const
{
  pnums.SetSize (NG_ELEMENT_MAXPOINTS);
  NG_ELEMENT_TYPE typ = Ng_GetElement (elnr+1, &pnums[0]);
  
  switch (typ)
    {
    case NG_TRIG: pnums.SetSize(3); break;
    case NG_TRIG6: pnums.SetSize(6); break;
    case NG_QUAD: pnums.SetSize(4); break;
    case NG_QUAD6: pnums.SetSize(6); break;

    case NG_TET: pnums.SetSize(4); break;
    case NG_TET10: pnums.SetSize(10); break;
    case NG_PYRAMID: pnums.SetSize(5); break;
    case NG_PRISM: pnums.SetSize(6); break;
    case NG_PRISM12: pnums.SetSize(12); break;
    case NG_HEX: pnums.SetSize(8); break;

    default:
      {
	stringstream err;
	err << "MeshAccess::GetElPNums undefined type "
	     << int(typ) << endl;
	throw Exception (err.str());
      }
    }

  for (int i = 0; i < pnums.Size(); i++)
    pnums[i]--;
}



void MeshAccess :: GetElVertices (int elnr, ARRAY<int> & vnums) const
{
  vnums.SetSize (NG_ELEMENT_MAXPOINTS);
  NG_ELEMENT_TYPE typ = Ng_GetElement (elnr+1, &vnums[0]);
  
  switch (typ)
    {
    case NG_TRIG: vnums.SetSize(3); break;
    case NG_TRIG6: vnums.SetSize(3); break;
    case NG_QUAD: vnums.SetSize(4); break;
    case NG_QUAD6: vnums.SetSize(4); break;

    case NG_TET: vnums.SetSize(4); break;
    case NG_TET10: vnums.SetSize(4); break;
    case NG_PYRAMID: vnums.SetSize(5); break;
    case NG_PRISM: vnums.SetSize(6); break;
    case NG_PRISM12: vnums.SetSize(6); break;
    case NG_HEX: vnums.SetSize(8); break;

    default:
      {
	stringstream err;
	err << "MeshAccess::GetElPNums undefined type "
	     << int(typ) << endl;
	throw Exception (err.str());
      }
    }

  for (int i = 0; i < vnums.Size(); i++)
    vnums[i]--;
}

  
void MeshAccess :: 
GetElEdges (int elnr, ARRAY<int> & ednums, ARRAY<int> & orient) const
{
  ednums.SetSize (12);
  orient.SetSize (12);
  int ned = 
    Ng_GetElement_Edges (elnr+1, &ednums[0], &orient[0]);
  ednums.SetSize (ned);
  orient.SetSize (ned);
  for (int i = 0; i < ned; i++)
    ednums[i]--;
}

void MeshAccess :: 
GetElEdges (int elnr, ARRAY<int> & ednums) const
{
  ednums.SetSize (12);
  int ned = 
    Ng_GetElement_Edges (elnr+1, &ednums[0], 0);
  ednums.SetSize (ned);
  for (int i = 0; i < ned; i++)
    ednums[i]--;
}

  
void MeshAccess :: 
GetSElEdges (int selnr, ARRAY<int> & ednums) const
{
  ednums.SetSize (4);
  int ned = 
    Ng_GetSurfaceElement_Edges (selnr+1, &ednums[0], 0);
  ednums.SetSize (ned);
  for (int i = 0; i < ned; i++)
    ednums[i]--;
}

void MeshAccess :: 
GetSElEdges (int selnr, ARRAY<int> & ednums, ARRAY<int> & orient) const
{
  ednums.SetSize (4);
  orient.SetSize (4);
  int ned = 
    Ng_GetSurfaceElement_Edges (selnr+1, &ednums[0], &orient[0]);
  ednums.SetSize (ned);
  orient.SetSize (ned);

  for (int i = 0; i < ned; i++)
    ednums[i]--;
}

  
void MeshAccess :: 
GetElFaces (int elnr, ARRAY<int> & fnums) const
{
  fnums.SetSize (6);
  int nfa = 
    Ng_GetElement_Faces (elnr+1, &fnums[0], 0);
  fnums.SetSize (nfa);

  for (int i = 0; i < nfa; i++)
    fnums[i]--;
}

void MeshAccess :: 
GetElFaces (int elnr, ARRAY<int> & fnums, ARRAY<int> & orient) const
{
  fnums.SetSize (6);
  orient.SetSize (6);
  int nfa = 
    Ng_GetElement_Faces (elnr+1, &fnums[0], &orient[0]);
  fnums.SetSize (nfa);
  orient.SetSize (nfa);

  for (int i = 0; i < nfa; i++)
    fnums[i]--;
}


int MeshAccess :: 
GetSElFace (int selnr) const
{
  return Ng_GetSurfaceElement_Face (selnr+1, 0)-1;
}
  
void MeshAccess :: 
GetSElFace (int selnr, int & fnum, int & orient) const
{
  fnum = Ng_GetSurfaceElement_Face (selnr+1, &orient);
  fnum--;
}

  
void MeshAccess :: GetFacePNums (int fnr, ARRAY<int> & pnums) const
{
  pnums.SetSize(4);
  int nv = Ng_GetFace_Vertices (fnr+1, &pnums[0]);
  pnums.SetSize(nv);
  for (int i = 0; i < nv; i++) 
    pnums[i]--;
}

 
void MeshAccess :: GetFaceEdges (int fnr, ARRAY<int> & edges) const
{
  edges.SetSize(4);
  int ned = Ng_GetFace_Edges (fnr+1, &edges[0]);
  edges.SetSize(ned);
  for (int i = 0; i < ned; i++) 
    edges[i]--;
}
 
void MeshAccess :: GetEdgePNums (int fnr, int & pn1, int & pn2) const
{
  int v2[2];
  Ng_GetEdge_Vertices (fnr+1, v2);
  pn1 = v2[0]-1;
  pn2 = v2[1]-1;
}




  
void MeshAccess :: GetSElPNums (int selnr, ARRAY<int> & pnums) const
{
  pnums.SetSize (NG_SURFACE_ELEMENT_MAXPOINTS);
  NG_ELEMENT_TYPE typ = Ng_GetSurfaceElement (selnr+1, &pnums[0]);
  switch (typ)
    {
    case NG_SEGM: pnums.SetSize(2); break;
    case NG_SEGM3: pnums.SetSize(3); break;

    case NG_TRIG: pnums.SetSize(3); break;
    case NG_TRIG6: pnums.SetSize(6); break;
    case NG_QUAD: pnums.SetSize(4); break;
    case NG_QUAD6: pnums.SetSize(6); break;

    default:
      {
	cerr << "MeshAccess::GetSElPNums undefined type "
	     << int(Ng_GetSurfaceElement (selnr, &pnums[0])) << endl;
      }

    }
  for (int i = 0; i < pnums.Size(); i++) pnums[i]--;
}


void MeshAccess :: GetSElVertices (int selnr, ARRAY<int> & vnums) const
{
  vnums.SetSize (NG_SURFACE_ELEMENT_MAXPOINTS);
  NG_ELEMENT_TYPE typ = Ng_GetSurfaceElement (selnr+1, &vnums[0]);
  switch (typ)
    {
    case NG_SEGM: vnums.SetSize(2); break;
    case NG_SEGM3: vnums.SetSize(2); break;

    case NG_TRIG: vnums.SetSize(3); break;
    case NG_TRIG6: vnums.SetSize(3); break;
    case NG_QUAD: vnums.SetSize(4); break;
    case NG_QUAD6: vnums.SetSize(4); break;

    default:
      {
	cerr << "MeshAccess::GetSElPNums undefined type "
	     << int(Ng_GetSurfaceElement (selnr+1, &vnums[0])) << endl;
      }

    }
  for (int i = 0; i < vnums.Size(); i++) vnums[i]--;
}


void MeshAccess :: GetVertexElements (int vnr, ARRAY<int> & elnrs) const
{
  int nel = Ng_GetNVertexElements (vnr+1);
  elnrs.SetSize (nel);
  Ng_GetVertexElements (vnr+1, &elnrs[0]);
  for (int j = 0; j < nel; j++)
    elnrs[j]--;
}






  


#ifdef NETGEN_ELTRANS

void MeshAccess ::
GetElementTransformation (int elnr, ElementTransformation & eltrans,
			  LocalHeap & lh) const
{
  int elind = Ng_GetElementIndex (elnr+1)-1;
  eltrans.SetElement (0, elnr, elind);
}

#else

void MeshAccess :: 
GetElementTransformation (int elnr, ElementTransformation & eltrans,
			  LocalHeap & lh) const
{
  int elind = Ng_GetElementIndex (elnr+1)-1;
  int pnums[NG_ELEMENT_MAXPOINTS];
  int np;

  NG_ELEMENT_TYPE et = Ng_GetElement (elnr+1, pnums);

  switch (et)
    {
    case NG_TRIG:
        eltrans.SetElement (&trig1, elnr, elind); break;
    case NG_TRIG6:
      eltrans.SetElement (&trig2, elnr, elind); break;
    case NG_QUAD:
      eltrans.SetElement (&quad1, elnr, elind); break;
    case NG_QUAD6:
      eltrans.SetElement (&quad2, elnr, elind); break;

    case NG_TET:
      eltrans.SetElement (&tet1, elnr, elind); break;
    case NG_TET10:
      eltrans.SetElement (&tet2, elnr, elind); break;
    case NG_PRISM:
      eltrans.SetElement (&prism1, elnr, elind); break;
    case NG_PRISM12:
      eltrans.SetElement (&prism2, elnr, elind); break;
    case NG_PYRAMID:
      eltrans.SetElement (&pyramid1, elnr, elind); break;
    case NG_HEX:
      eltrans.SetElement (&hex1, elnr, elind); break;

    default:
      cerr << "element transformation for element " << int(et)
	   << " not defined" << endl;
    }

  np = eltrans.GetElement().GetNDof();

  double point[3];
  int i, j;

  int dim = GetDimension();
  
  FlatMatrix<> & pmat = const_cast<FlatMatrix<>&> (eltrans.PointMatrix());
  pmat.AssignMemory (dim, np, lh);

  for (i = 0; i < np; i++)
    {
      Ng_GetPoint (pnums[i], point);
      for (j = 0; j < dim; j++)
	pmat(j, i) = point[j];
    }
}

#endif



#ifdef NETGEN_ELTRANS

void MeshAccess ::
GetSurfaceElementTransformation (int elnr, ElementTransformation & eltrans,
				 LocalHeap & lh) const
{
  int elind = Ng_GetSurfaceElementIndex (elnr+1)-1;
  eltrans.SetElement (1, elnr, elind);
}    

#else

void MeshAccess :: 
GetSurfaceElementTransformation (int elnr, ElementTransformation & eltrans, 
				 LocalHeap & lh) const
{
  try
    {
      int elind = Ng_GetSurfaceElementIndex (elnr+1)-1;
      
      int np;
      int pnums[NG_SURFACE_ELEMENT_MAXPOINTS];
      
      NG_ELEMENT_TYPE et = Ng_GetSurfaceElement (elnr+1, pnums);
      switch (et)
	{
	case NG_SEGM:
	  eltrans.SetElement (&segm1, elnr, elind); break;
	case NG_SEGM3:
	  eltrans.SetElement (&segm2, elnr, elind); break;
	  
	case NG_TRIG:
	  eltrans.SetElement (&trig1, elnr, elind); break;
	case NG_TRIG6:
	  eltrans.SetElement (&trig2, elnr, elind); break;
	case NG_QUAD:
	  eltrans.SetElement (&quad1, elnr, elind); break;
	case NG_QUAD6:
	  eltrans.SetElement (&quad2, elnr, elind); break;
	  
	default:
	  cerr << "surface element transformation for element " 
	       << int(et)
	       << " not defined" << endl;
	}

      np = eltrans.GetElement().GetNDof();

	
      FlatMatrix<> & pmat = const_cast<FlatMatrix<>&> (eltrans.PointMatrix());
      FlatMatrix<> & nvmat = const_cast<FlatMatrix<>&> (eltrans.NVMatrix());
      
      int dim = GetDimension();
      
      pmat.AssignMemory (dim, np, lh);
      nvmat.AssignMemory (dim, np, lh);
      
      double point[3], nv[3];
      int i, j;
      
      for (i = 0; i < np; i++)
	{
	  Ng_GetPoint (pnums[i], point);
	  for (j = 0; j < dim; j++)
	    pmat(j, i) = point[j];
	  
	  Ng_GetNormalVector (elnr+1, i+1, nv);
	  for (j = 0; j < dim; j++)
	    nvmat(j, i) = nv[j];
	}
      
      if (dim == 2)
	{
	  double tx = pmat(0,1)-pmat(0,0);
	  double ty = pmat(1,1)-pmat(1,0);
	  double len = sqrt (tx*tx+ty*ty);
	  if (len != 0)
	    { tx /= len; ty /= len; }
	  for (i = 0; i < np; i++)
	    {
	      nvmat(0,i) = ty;
	      nvmat(1,i) = -tx;
	    }
	}
    }
  catch (Exception & e)
    {
      stringstream ost;
      ost << "in GetSurfaceElementTransformation" << endl;
      e.Append (ost.str());
      throw;
    }
  catch (exception & e)
    {
      throw (Exception (string(e.what()) +
			string("\n in GetSurfaceElementTransformation\n")));
    }

}

#endif



  
double MeshAccess :: ElementVolume (int elnr) const
{

  const FiniteElement * fe;
  switch (GetElType (elnr))
    {
    case ET_TRIG: fe = &trig0; break;
    case ET_QUAD: fe = &quad0; break;
    case ET_TET: fe = &tet0; break;
    case ET_PYRAMID: fe = &pyramid0; break;
    case ET_PRISM: fe = &prism0; break;
    case ET_HEX: fe = &hex0; break;
    default:
      {
	cerr << "ElementVolume not implemented for el " << GetElType(elnr) << endl;
      }
    }
  
  char d[10000];
  LocalHeap lh(d, 10000);

  ElementTransformation trans;
  GetElementTransformation (elnr, trans, lh);
  ConstantCoefficientFunction ccf(1);


  if (GetDimension() == 2)
    {
      SourceIntegrator<2> si( &ccf );
      FlatVector<> elvec;
      si.AssembleElementVector (*fe, trans, elvec, lh);
      return elvec(0);
    }
  else
    {
      SourceIntegrator<3> si( &ccf );
      FlatVector<> elvec;
      si.AssembleElementVector (*fe, trans, elvec, lh);
      return elvec(0);
    }
}


  
double MeshAccess :: SurfaceElementVolume (int selnr) const
{
  const FiniteElement * fe;
  switch (GetSElType (selnr))
    {
    case ET_TRIG: fe = &trig0; break;
    case ET_QUAD: fe = &quad0; break;
    default:
      {
	cerr << "SurfaceElementVolume not implemented for el " << GetElType(selnr) << endl;
	return 0;
      }
    }

  char d[10000];
  LocalHeap lh(d, 10000);

  ElementTransformation trans;
  GetSurfaceElementTransformation (selnr, trans, lh);
  ConstantCoefficientFunction ccf(1);

  if (GetDimension() == 2)
    {
      NeumannIntegrator<2> si( &ccf );
      FlatVector<> elvec;
      si.AssembleElementVector (*fe, trans, elvec, lh);
      return elvec(0);
    }
  else
    {
      NeumannIntegrator<3> si( &ccf );
      FlatVector<> elvec;
      si.AssembleElementVector (*fe, trans, elvec, lh);
      return elvec(0);
    }
}

int MeshAccess :: GetNLevels() const
{
  return Ng_GetNLevels();
}


int MeshAccess :: FindElementOfPoint (FlatVector<double> point,
				      IntegrationPoint & ip,
				      bool build_searchtree,
				      int index) const
{
  int elnr;
  double lami[3];
  elnr = Ng_FindElementOfPoint (&point(0), lami, build_searchtree, index);

  if (elnr == 0) return -1;

  /*
  if (fel.ElementType() == ET_TRIG)
    {
      double hlam1 = lami[0];
      double hlam2 = lami[1];
      lami[0] = 1-hlam1-hlam2;
      lami[1] = hlam1;
    }
  */

  // if (fel.ElementType() == ET_TET)
  if (0)
    {
      ip(0) = 1-lami[0]-lami[1]-lami[2];
      ip(1) = lami[0];
      ip(2) = lami[1];
      /*
      double hlam1 = lami[0];
      double hlam2 = lami[1];
      double hlam3 = lami[2];
      lami[0] = 1-hlam1-hlam2-hlam3;
      lami[1] = hlam1;
      lami[2] = hlam2;
      */
    }
  else
    for (int k = 0; k < 3; k++)
      ip(k) = lami[k];

  return elnr-1;
}



int MeshAccess :: GetNPairsPeriodicVertices () const 
{
  return Ng_GetNPeriodicVertices();
}
 
void MeshAccess :: GetPeriodicVertices (ARRAY<INT<2> > & pairs) const
{
  int npairs;

  npairs = Ng_GetNPeriodicVertices ();
  pairs.SetSize (npairs);

  Ng_GetPeriodicVertices (&pairs[0][0]);
  for (int i = 0; i < pairs.Size(); i++)
    {
      pairs[i][0]--;
      pairs[i][1]--;
    }
}

int MeshAccess :: GetNPairsPeriodicEdges () const 
{
  return Ng_GetNPeriodicEdges();
}
 
void MeshAccess :: GetPeriodicEdges (ARRAY<INT<2> > & pairs) const
{
  int npairs;

  npairs = Ng_GetNPeriodicEdges ();
  pairs.SetSize (npairs);

  Ng_GetPeriodicEdges (&pairs[0][0]);
  for (int i = 0; i < pairs.Size(); i++)
    {
      pairs[i][0]--;
      pairs[i][1]--;
    }
}

/*
///// Added by Roman Stainko ....
void MeshAccess::GetVertexElements( int vnr, ARRAY<int>& elems) const
{
  int nel = Ng_GetVertex_NElements( vnr+1 );
  elems.SetSize( nel );
  
  Ng_GetVertex_Elements( vnr+1, &elems[0] );
  
  for( int i=0; i<nel; i++ )
    elems[i]--;
 
}
*/

///// Added by Roman Stainko ....
void MeshAccess::GetVertexSurfaceElements( int vnr, ARRAY<int>& elems) const
{
  int nel = Ng_GetVertex_NSurfaceElements( vnr+1 );
  elems.SetSize( nel );
  
  Ng_GetVertex_SurfaceElements( vnr+1, &elems[0] );

  for( int i=0; i<nel; i++ )
    elems[i]--;
}
}