#ifndef __MESH_H
#define __MESH_H
/*!
  \file    mesh.hh
  \brief   2D triangular mesh with segmented boundary
  \date    30 April 2010
  \author  Clemens Pechstein

  Institute of Computational Mathematics
  Johannes Kepler University Linz, Austria

  Tutorials -- Numerical methods for elliptic PDEs
  summer semester 2010
*/
#include <vector>
#include <valarray>
#include <iostream>
#include <string>
#include "vec.hh"
#include "sparsematrix.hh"

using namespace std;

//from vec.hh:
//typedef Vec<2, double> Point2D;

/// triangle: 3 node indices
typedef Vec<3, int> Trig;

/// segment: 2 node indices
typedef Vec<2, int> Segment;

/// boundary condition: 
enum BoundaryCondition {
  BC_INTERNAL  = 0, ///< none
  BC_DIRICHLET = 1, 
  BC_NEUMANN   = 2,
  BC_ROBIN     = 3
};

inline
bool isDirichletOrRobin (BoundaryCondition bc)
{
  return (bc == BC_DIRICHLET) || (bc == BC_ROBIN); 
}

inline
string BC_STR(BoundaryCondition bc)
{
  if (bc == BC_INTERNAL) return "internal";
  else if (bc == BC_DIRICHLET) return "Dirichlet";
  else if (bc == BC_NEUMANN) return "Neumann";
  else if (bc == BC_ROBIN) return "Robin";
  else return "unknown";
}



/// triangular 2D mesh with segmented boundary
class Mesh
{
public:
  //--- public members 
  
  int nVerts;    ///< number of vertices
  int nTrigs;    ///< number of triangles
  int nSegments; ///< number of segments on the boundary
  
  vector<Point2D> verts;    ///< coordinates of the vertices
  vector<Trig> trigs;       ///< triangles (3 vertex indices)
  vector<Segment> segments; ///< boundary segments (2 vertex indices)
  
  vector<BoundaryCondition> bcSegments;  ///< boundary conditions at the segments
  vector<BoundaryCondition> bcVerts;     ///< boundary conditions at (all) the vertices

  //--- public functions
  
  Mesh () : nVerts(0), nTrigs(0), nSegments(0),
	    verts(), trigs(), segments(), bcSegments(), bcVerts()
  { ; }

  ~Mesh () { ; }
  
  /// reserve memory: for  nV vertices,  nT triangles,  nS segments
  void reserve (int nV, int nT, int nS=0);
  
  /// set boundary condition on segment and auto-set boundary conditions on the two vertices
  void setSegmentBC (int s, BoundaryCondition bc);
  
  /// get FEM-matrix shape for this mesh
  void getMatrixShape (SparseShape& ss) const;
  
  /// create refined mesh out of this one
  void refineUniform (Mesh& mesh) const;
  
  /// create refined mesh
  inline
  void refineUniform ()
  { Mesh tmp; tmp.copyFrom (*this); tmp.refineUniform (*this); }
  
  /**
   * matlab output
   * write mesh (and solution data) and two plot commands into a matlab-file
   */
  void matlabOutput (const char* filename, const Vector& solution=Vector()) const;
  
  /**
   * write mesh (+ solution) to a .gmv-file
   * GMV = General Mesh Viewer, see
   *
   * http://www-xdiv.lanl.gov/XCM/gmv/GMVHome.html
   *
   * @param  filename  name of the .gmv-file
   * @param  solution  optional: field of size nVerts
   * @param  z         optional: if true plot field above as z-axis
   */
  void gmvOutput (const char* filename, const Vector& solution=Vector(), bool z=false,
                  const Vector& error=Vector()) const;

  /// gnuplot output
  void gnuplotOutput (const char* fileprefix, const Vector& solution) const;

  /// output information
  void report (ostream& os) const;

  /// get longest element edge = largest triangle diameter
  double longestEdge () const;
  
  /// copy from other mesh
  void copyFrom (const Mesh& m);

private:
  
  Mesh(const Mesh& mesh);            ///< hidden in order to prevent unconcious copy: use "copyFrom" to copy !
  Mesh& operator=(const Mesh& mesh); ///< hidden in order to prevent unconcious copy: use "copyFrom" to copy !
  
}; // class Mesh


inline
void refineUniform (const Mesh& coarse, Mesh& fine, int steps=1)
{
  if (steps == 0)
    fine.copyFrom (coarse);
  if (steps == 1)
    coarse.refineUniform (fine);
  else if (steps > 1)
    {
      Mesh medium;
      coarse.refineUniform (medium);
      refineUniform (medium, fine, steps-1);
    }
}


#endif // __MESH_H