/*!
  \file    meshdemo.cc
  \brief   demo for mesh
  \date    4 May 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"

using namespace std;

int main ()
{
  Mesh mesh;
  
  mesh.reserve (9, 8, 8); // 9 nodes(vertices), 8 elements, and 8 boundary segments
  
  // set vertex coordinates
  mesh.verts[0].set (0.0, 0.0);
  mesh.verts[1].set (0.5, 0.0);
  mesh.verts[2].set (1.0, 0.0);
  mesh.verts[3].set (0.0, 0.5);
  mesh.verts[4].set (0.5, 0.5);
  mesh.verts[5].set (1.0, 0.5);
  mesh.verts[6].set (0.0, 1.0);
  mesh.verts[7].set (0.5, 1.0);
  mesh.verts[8].set (1.0, 1.0);
  
  // create triangles by 3 vertex indices each 
  mesh.trigs[0].set (0, 4, 3);   // vertex 0: (0.0, 0.0), 4: (0.5, 0.5), and 3: (0.0, 0.5)
  mesh.trigs[1].set (0, 1, 4);   // vertex 0: (0.0, 0.0), 1: (0.5, 0.0), and 4: (0.5, 0.5)
  mesh.trigs[2].set (1, 5, 4);
  mesh.trigs[3].set (1, 2, 5);
  mesh.trigs[4].set (3, 7, 6);
  mesh.trigs[5].set (3, 4, 7);
  mesh.trigs[6].set (4, 8, 7);
  mesh.trigs[7].set (4, 5, 8);
  
  // create boundary segments by 2 vertex indices each
  mesh.segments[0].set (0, 1); mesh.setSegmentBC (0, BC_NEUMANN);
  mesh.segments[1].set (1, 2); mesh.setSegmentBC (1, BC_NEUMANN);
  mesh.segments[2].set (2, 5); mesh.setSegmentBC (2, BC_NEUMANN);
  mesh.segments[3].set (5, 8); mesh.setSegmentBC (3, BC_NEUMANN);
  mesh.segments[4].set (8, 7); mesh.setSegmentBC (4, BC_NEUMANN);
  mesh.segments[5].set (7, 6); mesh.setSegmentBC (5, BC_NEUMANN);
  mesh.segments[6].set (6, 3); mesh.setSegmentBC (6, BC_NEUMANN);
  mesh.segments[7].set (3, 0); mesh.setSegmentBC (7, BC_NEUMANN);
  
  cout << "mesh info:" << endl
       << "==========" << endl;

  mesh.report (cout);
  cout << "h = " << mesh.longestEdge() << endl << endl;

  cout << "coordinates of vertex 4: "
       << mesh.verts[4] << endl << endl;

  cout << "vertex indices of triangle 2: "
       << mesh.trigs[2] << endl << endl;

  cout << "coordinates of triangle 2: "
       << mesh.verts[mesh.trigs[2][0]] << ", "
       << mesh.verts[mesh.trigs[2][1]] << ", "
       << mesh.verts[mesh.trigs[2][2]] << endl << endl;
  
  mesh.refineUniform (); // perform one step of uniform refinement
  
  cout << "mesh after refinement:" << endl
       << "======================" << endl;
  mesh.report (cout);
  cout << "h = " << mesh.longestEdge() << endl << endl;
  
  cout << "writing mesh to 'output.m'..." << endl;
  mesh.matlabOutput ("output.m");
  
  // CLEMENS
  cout << "writing mesh to 'output.gmv'..." << endl;
  mesh.gmvOutput ("output.gmv");
  
  return 0;
}