//********************************************************************************
//*
//*  C++ finite element method for heat equation
//*  James sandham
//*  15 April 2015
//*
//********************************************************************************

//********************************************************************************
//
// HeatFE is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
//********************************************************************************

#include<iostream>
#include<fstream>
#include<stdlib.h>
#include<cstring>
#include"Mesh.h"
using std::cout;
using std::endl;
using std::ifstream;

const int MAX_NUM_GROUPS = 20;
const int MAX_NUM_ELEM_TYP = 20;
const int MAX_CHARS_PER_LINE = 512;
const int MAX_TOKENS_PER_LINE = 20;
const char* DELIMITER = " ";



//----------------------------------------------------------------------------
// Constructor for Mesh
//----------------------------------------------------------------------------
Mesh::Mesh(const char *file)
{
  //create a file-reading object
  ifstream fin;
  fin.open(file);

  int flag=0,index=0,typ=0,grp=0;
  int tl[MAX_NUM_ELEM_TYP]={}; //list of element types
  int gl[MAX_NUM_GROUPS]={};   //list of groups

  Dim = 0;     //dimension of problem
  Ng = 0;      //number of element groups
  N = 0;       //number of points 
  Nte = 0;     //total number of elements 
  Ne = 0;      //number of interior elements
  Ne_b = 0;    //number of boundary elements
  Npe = 0;     //number of points per interior element
  Npe_b = 0;   //number of points per boundary element
  Type = 0;    //interior element type
  Type_b = 0;  //boundary element type

  //scan through file (first pass)
  while(!fin.eof())
  {
    char buf[MAX_CHARS_PER_LINE];
    fin.getline(buf,MAX_CHARS_PER_LINE);

    const char* token[MAX_TOKENS_PER_LINE]={};

    int n = 0;

    // parse the line
    token[0] = strtok(buf,DELIMITER);
    if(token[0]){
      if(strcmp(token[0],"$Nodes")==0) {flag=1; index=-2;}
      if(strcmp(token[0],"$Elements")==0) {flag=2; index=-2;}
      if(strcmp(token[0],"$EndNodes")==0) {flag=0; index=-2;}
      if(strcmp(token[0],"$EndElements")==0) {flag=0; index=-2;}
      for(n=1;n<MAX_TOKENS_PER_LINE;n++){
        token[n] = strtok(0,DELIMITER);
        if(!token[n]) break;
      }
    }

    //process the line
    if(index==-1){
      if(flag==1) {N = atoi(token[0]);}
      if(flag==2) {Nte = atoi(token[0]);}
    }
    else if(index>-1){
      if(flag==2){
        typ = atoi(token[1]);  //element type
        grp = atoi(token[3]);  //group
        if(typ==15){           //15 corresponds to a point
        }
        else if(typ>Type){
          Type=typ;
          Ne = 1;
        }
        else if(typ==Type){
          Ne++;
        }
        for(int i=0;i<MAX_NUM_GROUPS;i++){
          if(gl[i]==0) {gl[i]=grp; Ng++; break;}
          if(gl[i]==grp) {break;}
        }
        for(int i=0;i<MAX_NUM_ELEM_TYP;i++){
          if(tl[i]==0) {tl[i]=typ; break;}
          if(tl[i]==typ) {break;}
        }
      }
    }
    index++;
  }
  Ne_b = Nte-Ne;

  //initialize model
  groups = new int[Ng];
  xpoints = new double[N];
  ypoints = new double[N];
  zpoints = new double[N];
  connect = new int*[Ne];
  connect_b = new int*[Ne_b];

  for(int i=0;i<Ng;i++){groups[i] = gl[i];}

  switch (Type)
  {
    case 1:      //linear 1D lines
      Npe = 2;
      Npe_b = 1;
      Type_b = 15;
      Dim = 1;
      break;
    case 2:      //linear 2D triangles
      Npe = 3;
      Npe_b = 2;
      Type_b = 1;
      Dim = 2;
      break;
    case 3:      //linear 2D quadrangles
      Npe = 4;
      Npe_b = 2;
      Type_b = 1;
      Dim = 2;
      break;
    case 4:      //linear 3D tetrahedra
      Npe = 4;
      Npe_b = 3;
      Type_b = 2;
      Dim = 3;
      break;
    case 8:      //quadratic 1D lines
      Npe = 3;
      Npe_b = 1;
      Type_b =15;
      Dim = 1;
      break;
    case 9:      //quadratic 2D triangles
      Npe = 6;
      Npe_b = 3;
      Type_b = 8;
      Dim = 2;
      break;
    case 10:     //quadratic 2D quadrangles
      Npe = 8;
      Npe_b = 3;
      Type_b = 8;
      Dim = 2;
      break;
    case 11:     //quadratic 3D tetrahedra
      Npe = 10;
      Npe_b = 6;
      Type_b = 9;
      Dim = 3;
      break;
  }
  for(int i=0;i<Ne;i++)
    connect[i] = new int[Npe];
  for(int i=0;i<Ne_b;i++)
    connect_b[i] = new int[Npe_b+1];  //include extra column for group numbers


  //return to beginning of file
  fin.clear();
  fin.seekg(0,fin.beg);

  //scan through file (second pass)
  index=0; flag=0;
  while(!fin.eof())
  {
    char buf[MAX_CHARS_PER_LINE];
    fin.getline(buf,MAX_CHARS_PER_LINE);

    const char* token[MAX_TOKENS_PER_LINE]={};

    int n = 0;

    // parse the line
    token[0] = strtok(buf,DELIMITER);
    if(token[0]){
      if(strcmp(token[0],"$Nodes")==0) {flag=1; index=-2;}
      if(strcmp(token[0],"$Elements")==0) {flag=2; index=-2;}
      if(strcmp(token[0],"$EndNodes")==0) {flag=0; index=-2;}
      if(strcmp(token[0],"$EndElements")==0) {flag=0; index=-2;}
      for(n=1;n<MAX_TOKENS_PER_LINE;n++){
        token[n] = strtok(0,DELIMITER);
        if(!token[n]) break;
      }
    }

    //process the line (fill model arrays)
    if(index>-1){
      if(flag==1){
        xpoints[index] = strtod(token[1],NULL);
        ypoints[index] = strtod(token[2],NULL);
        zpoints[index] = strtod(token[3],NULL);
      }
      if(flag==2){
        if(atoi(token[1])==Type){
          switch (Type)
          {
          case 1:       //2-point 1D line
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            break;
          case 2:       //3-point 2D triangle
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            break;
          case 3:       //4-point 2D quadrangle
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            connect[index-Ne_b][3] = atoi(token[8]);
            break;
          case 4:       //4-point 3D tetrahedra
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            connect[index-Ne_b][3] = atoi(token[8]);
            break;
          case 8:       //3-point 1D line
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            break;
          case 9:       //6-point 2D triangle
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            connect[index-Ne_b][3] = atoi(token[8]);
            connect[index-Ne_b][4] = atoi(token[9]);
            connect[index-Ne_b][5] = atoi(token[10]);
            break;
          case 10:      //8-point 2D quadrangle
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            connect[index-Ne_b][3] = atoi(token[8]);
            connect[index-Ne_b][4] = atoi(token[9]);
            connect[index-Ne_b][5] = atoi(token[10]);
            connect[index-Ne_b][6] = atoi(token[11]);
            connect[index-Ne_b][7] = atoi(token[12]);
            break;
          case 11:      //10-point 3D tetrahedra
            connect[index-Ne_b][0] = atoi(token[5]);
            connect[index-Ne_b][1] = atoi(token[6]);
            connect[index-Ne_b][2] = atoi(token[7]);
            connect[index-Ne_b][3] = atoi(token[8]);
            connect[index-Ne_b][4] = atoi(token[9]);
            connect[index-Ne_b][5] = atoi(token[10]);
            connect[index-Ne_b][6] = atoi(token[11]);
            connect[index-Ne_b][7] = atoi(token[12]);
            connect[index-Ne_b][8] = atoi(token[13]);
            connect[index-Ne_b][9] = atoi(token[14]);
            break;
          }
        }
        else{
          switch (Type)
          {
          case 1:       //1-point (2-point 1D line)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            break;
          case 2:       //2-point line (3-point 2D triangle)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            connect_b[index][2] = atoi(token[6]);
            break;
          case 3:       //2-point line (4-point 2D quadrangle)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            connect_b[index][2] = atoi(token[6]);
            break;
          case 4:       //3-point triangle (4-point 3D tetrahedra)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            connect_b[index][2] = atoi(token[6]);
            connect_b[index][3] = atoi(token[7]);
            break;
          case 8:       //1-point (3-point 1D line)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            break;
          case 9:       //3-point line (6-point 2D triangle)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            connect_b[index][2] = atoi(token[6]);
            connect_b[index][3] = atoi(token[7]);
            break;
          case 10:      //3-point line (8-point 2D quadrangle)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            connect_b[index][2] = atoi(token[6]);
            connect_b[index][3] = atoi(token[7]);
            break;
          case 11:      //6-point triangle (10-point 3D tetrahedra)
            connect_b[index][0] = atoi(token[3]);
            connect_b[index][1] = atoi(token[5]);
            connect_b[index][2] = atoi(token[6]);
            connect_b[index][3] = atoi(token[7]);
            connect_b[index][4] = atoi(token[8]);
            connect_b[index][5] = atoi(token[9]);
            connect_b[index][6] = atoi(token[10]);
            break;
          }
        }
      }
    }
    index++;
  }
  fin.close();
}



//----------------------------------------------------------------------------
// Delete Mesh
//----------------------------------------------------------------------------
void Mesh::deleteMesh()
{
  delete [] xpoints;
  delete [] ypoints;
  delete [] zpoints;

  for(int i=0;i<Ne;i++){
    delete [] connect[i];
  }
  delete [] connect;

  for(int i=0;i<Ne_b;i++){
    delete [] connect_b[i];
  }
  delete [] connect_b;
}
