/*!
  \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 "mesh.hh"
#include <fstream>
#include <string>
#include <map>

using namespace std;

void Mesh :: reserve (int nV, int nT, int nS)
{
  nVerts = nV;
  nTrigs = nT;
  nSegments = nS;
    
  verts.resize (nVerts);
  trigs.resize (nTrigs);
  segments.resize (nSegments);
  bcSegments.resize (nSegments);
  bcVerts.resize (nVerts);
  for (int i=0; i<nVerts; ++i) bcVerts[i] = BC_INTERNAL;
  for (int i=0; i<nSegments; ++i) bcSegments[i] = BC_INTERNAL;
}



void Mesh :: setSegmentBC (int s, BoundaryCondition bc)
{
  bcSegments[s] = bc;
  if (! isDirichletOrRobin (bcVerts[segments[s][0]])) 
    bcVerts[segments[s][0]] = bc;
  if (! isDirichletOrRobin (bcVerts[segments[s][1]])) 
    bcVerts[segments[s][1]] = bc;
}



void Mesh :: getMatrixShape (SparseShape& ss) const
{
  ss.resize (nVerts, nVerts);
  
  // PUT YOUR CODE HERE
}


/// ordered pair of two ints
struct Edge
{
  int n[2];
  
  Edge(int n1=0, int n2=0) { n[0] = min (n1, n2); n[1] = max (n1, n2); }
  Edge(const Edge& e) { n[0] = e.n[0]; n[1] = e.n[1]; }
  Edge& operator= (const Edge& e) { n[0] = e.n[0]; n[1] = e.n[1]; return *this;}
  
  const int operator[] (int i) const { assert (0 <= i && i < 2); return n[i]; }
  int& operator[] (int i) { assert (0 <= i && i < 2); return n[i]; }

  void set (int n1, int n2) { n[0] = min (n1, n2); n[1] = max (n1, n2); }

  bool operator== (const Edge& e) const { return (n[0] == e.n[0]) && (n[1] == e.n[1]); }
  bool operator< (const Edge& e) const { return (n[0] < e.n[0] || (n[0] == e.n[0] && n[1] < e.n[1])); }

}; // struct Edge



void Mesh :: refineUniform (Mesh& mesh) const
{
  typedef Edge KeyType;
  typedef map<KeyType, int> MapType;
  typedef pair<KeyType, int> PairType;
  MapType edgeInfo; // map edge -> new vertex index
  MapType::iterator it;
  int count; 
  PairType p;
  KeyType key;
  
  // generate the set of edges by inserting (duplicates won't be inserted twice!!!)
  for (int t=0; t<nTrigs; ++t)
    {
      p.first.set (trigs[t][0], trigs[t][1]);
      edgeInfo.insert (p);
      p.first.set (trigs[t][1], trigs[t][2]);
      edgeInfo.insert (p);
      p.first.set (trigs[t][2], trigs[t][0]);
      edgeInfo.insert (p);
    }
  
  // count the edges and at the same time give new vertex numbers
  count = nVerts;
  for (it=edgeInfo.begin (); it != edgeInfo.end (); ++it)
    {
      it->second = count;
      ++count;
    }
  
  // reserve for new mesh
  mesh.reserve (count, 4*nTrigs, 2*nSegments);
  
  // copy old verts
  for (int i=0; i<nVerts; ++i)
    {
      mesh.verts[i] = verts[i];
      mesh.bcVerts[i] = bcVerts[i];
    }
  // set new verts
  for (it = edgeInfo.begin (); it != edgeInfo.end (); ++it)
    {
      mesh.verts[it->second] = 0.5 * (verts[it->first[0]] + verts[it->first[1]]);
      mesh.bcVerts[it->second] = BC_INTERNAL;
    }

  // refine triangles
  int new01, new12, new20;
  count = 0;
  
  /*
       n[2]
        *
       / \ 
      /   \
 n20 x-----x n12
    / \   / \
   /   \ /   \
  *-----x-----*
 n[0]   n01    n[1]

   */

  for (int t=0; t<nTrigs; ++t)
    {
      key.set (trigs[t][0], trigs[t][1]);
      it = edgeInfo.find (key);
      new01 = it->second;
      
      key.set (trigs[t][1], trigs[t][2]);
      it = edgeInfo.find (key);
      new12 = it->second;

      key.set (trigs[t][2], trigs[t][0]);
      it = edgeInfo.find (key);
      new20 = it->second;
      
      mesh.trigs[count].set (trigs[t][0], new01, new20); ++count;
      mesh.trigs[count].set (trigs[t][1], new12, new01); ++count;
      mesh.trigs[count].set (trigs[t][2], new20, new12); ++count;
      mesh.trigs[count].set (new01, new12, new20); ++count;
    }
  
  count = 0;

  // refine segments
  for (int s=0; s<nSegments; ++s)
    {
      key.set (segments[s][0], segments[s][1]);
      it = edgeInfo.find (key);
      new01 = it->second;
      
      mesh.segments[count].set (segments[s][0], new01);
      mesh.bcSegments[count] = bcSegments[s];
      ++count;
      mesh.segments[count].set (new01, segments[s][1]);
      mesh.bcSegments[count] = bcSegments[s];
      ++count;
      mesh.bcVerts[new01] = bcSegments[s];
    }
  
} // Mesh :: refineUniform



void Mesh :: gmvOutput (const char* filename, const Vector& solution, bool z, const Vector& err) const
{
  bool ok = ((int)solution.size () == nVerts);

  ofstream ofs;
  ofs.open (filename);
  
  ofs << "gmvinput ascii" << endl << endl;
  
  ofs << "nodev " << nVerts << endl;
  for (int i=0; i<nVerts; ++i)
    {
      ofs << verts[i][0] << "\t" << verts[i][1] << "\t"
	  << ((ok && z) ? solution[i] : 0.0) << endl;
    }
  ofs << endl;
  
  ofs << "cells " << nTrigs + nSegments << endl;
  for (int t=0; t<nTrigs; ++t)
    ofs << "tri 3 " << trigs[t][0]+1 << " " << trigs[t][1]+1 << " " << trigs[t][2]+1 << endl;  
  for (int s=0; s<nSegments; ++s)
    ofs << "line 2 " << segments[s][0]+1 << " " << segments[s][1]+1 << endl;
  ofs << endl;
  
  ofs << "material 4 0" << endl; // 0..cell material
  ofs << "internal Dirichl Neumann Robin" << endl;
  for (int t=0; t<nTrigs; ++t) ofs << "1 ";
  ofs << endl;
  for (int s=0; s<nSegments; ++s)
    if (bcSegments[s] == BC_DIRICHLET) ofs << "2 ";
    else if (bcSegments[s] == BC_NEUMANN) ofs << "3 ";
    else if (bcSegments[s] == BC_ROBIN) ofs << "4 ";
    else ofs << "1 ";
  ofs << endl << endl;
  
  if (ok)
    {
      ofs << "variable" << endl;
      ofs << "solution 1" << endl; // 1..node field
      for (int i=0; i<nVerts; ++i) ofs << solution[i] << " ";
      ofs << endl;
      if (err.size() == nTrigs)
        {
          ofs << "error 0" << endl; //0..cell field
          for (int i=0; i<nTrigs; ++i) ofs << err[i] << " ";
          ofs << endl;
        }
      ofs << "endvars" << endl;
    }
  
  ofs << "endgmv" << endl;
  
  ofs.close ();
  
} // Mesh :: gmvOutput


void Mesh :: gnuplotOutput (const char* fileprefix, const Vector& solution) const
{
  string datafile = string(fileprefix) + ".gpd";
  ofstream ofs;
  ofs.open (datafile.c_str());
  for (int i=0; i<nVerts; ++i)
    {
      ofs << verts[i][0] << '\t' << verts[i][1] << '\t' << solution[i] << endl;
    }
  ofs.close ();
  
  datafile = string(fileprefix) + ".gp";
  ofs.open (datafile.c_str());
  ofs << "splot '" << fileprefix << ".gpd' using 1:2:3 with points" << endl
      << endl;
  ofs << "pause -1" << endl;
  ofs.close ();
  
} // Mesh :: gnuplotOutput



void Mesh :: matlabOutput (const char* filename, const Vector& solution) const
{
  ofstream ofs (filename);

  ofs << "p = [[" << verts[0][0] << "; " << verts[0][1] << "]";
  for (int i=1; i<nVerts; ++i)
    ofs << ", [" <<  verts[i][0] << "; " << verts[i][1] << "]";
  ofs << "];" << endl;
  ofs << endl;

  if (nSegments == 0)
    {
      ofs << "e = []" << endl;
    }
  else
    {
      ofs << "e = [[" << segments[0][0]+1 << "; " << segments[0][1]+1 << "]";
      for (int i=1; i<nSegments; ++i)
        ofs << ", [" << segments[i][0]+1 << "; " << segments[i][1]+1 << "]";
      ofs << "];" << endl;
      ofs << endl;
    }
  
  ofs << "t = [[" << trigs[0][0]+1 << "; " << trigs[0][1]+1 << "; " << trigs[0][2]+1 << "; 1]";
  for (int i=1; i<nTrigs; ++i)
    ofs << ", [" << trigs[i][0]+1 << "; " << trigs[i][1]+1 << "; " << trigs[i][2]+1 << "; 1]";
  ofs << "];" << endl;
  ofs << endl;
  
  if (solution.size() != nVerts)
    {
      ofs << "pdemesh (p, e, t);" << endl << endl;
    }
  else
    {
      ofs << "u = [" << solution[0];
      for (int i=1; i<nVerts; ++i)
        ofs << ", " << solution[i];
      ofs << "];" << endl;
      ofs << endl;
      
      ofs << "pdemesh (p, e, t, u);" << endl << endl;

      ofs << "pdesurf (p, t, transpose(u));" << endl << endl;
    }
  
} // Mesh :: matlabOutput



void Mesh :: report (ostream& os) const
{
  os << "VERTS [ " << nVerts << " ]" << endl;
  for (int i=0; i<nVerts; ++i)
    {
      os << i << ": " << verts[i] << "\t" << BC_STR(bcVerts[i]) << endl;
    }
  os << endl;
  
  os << "TRIGS [ " << nTrigs << " ]" << endl;
  for (int i=0; i<nTrigs; ++i)
    {
      os << i << ": " << trigs[i] << endl;
    }
  os << endl;
  
  os << "SEGMENTS [ " << nSegments << " ]" << endl;
  for (int i=0; i<nSegments; ++i)
    {
      os << i << ": " << segments[i] << "\t" << BC_STR(bcSegments[i]) << endl;
    }
  os << endl;

} // Mesh :: report



double Mesh :: longestEdge () const
{
  double longest = 0;
  Vec<3, Point2D> p;
  
  for (int t=0; t<nTrigs; ++t)
    {
      for (int j=0; j<3; ++j) p[j] = verts[trigs[t][j]];
      longest = max(longest, l2norm (p[1] - p[0]));
      longest = max(longest, l2norm (p[2] - p[0]));
      longest = max(longest, l2norm (p[2] - p[1]));
    }
  
  return longest;
  
} // Mesh :: longestEdge



void Mesh :: copyFrom (const Mesh& m)
{
  nVerts = m.nVerts;
  nTrigs = m.nTrigs;
  nSegments = m.nSegments;
  verts = m.verts;
  trigs = m.trigs;
  segments = m.segments;
  bcSegments = m.bcSegments;
  bcVerts = m.bcVerts;
  
} // Mesh :: copyFrom


// EOF