//********************************************************************************
//*
//*  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"FeModel.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 Finite Element Model (read mesh)
//-----------------------------------------------------------------------------
FeModel::FeModel(char* file)
{
  //create a file-reading object
  ifstream fin;
  fin.open(file);
  model = 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
  N = 0;     //number of points 
  Nte = 0;   //total number of elements 
  Ne = 0;    //number of interior elements
  Nb = 0;    //number of boundary elements
  Ng = 0;    //number of groups
  Npe = 0;   //number of points per interior element
  Nbpe = 0;  //number of points per boundary element
  Type = 0;  //interior element type
  bType = 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++;
  }
  Nb = Nte-Ne;

  //initialize model
  xpoints = new double[N];
  ypoints = new double[N];
  zpoints = new double[N];
  connect = new int*[Ne];
  bconnect = new int*[Nb];
  switch (Type)
  {
    case 1:      //linear 1D lines
      Npe = 2;
      Nbpe = 1;
      bType = 15;
      Dim = 1;
      break;
    case 2:      //linear 2D triangles
      Npe = 3;
      Nbpe = 2;
      bType = 1;
      Dim = 2;
      break;
    case 3:      //linear 2D quadrangles
      Npe = 4;
      Nbpe = 2;
      bType = 1;
      Dim = 2;
      break;
    case 4:      //linear 3D tetrahedra
      Npe = 4;
      Nbpe = 3;
      bType = 2;
      Dim = 3;
      break;
    case 8:      //quadratic 1D lines
      Npe = 3;
      Nbpe = 1;
      bType =15;
      Dim = 1;
      break;
    case 9:      //quadratic 2D triangles
      Npe = 6;
      Nbpe = 3;
      bType = 8;
      Dim = 2;
      break;
    case 10:     //quadratic 2D quadrangles
      Npe = 8;
      Nbpe = 3;
      bType = 8;
      Dim = 2;
      break;
    case 11:     //quadratic 3D tetrahedra
      Npe = 10;
      Nbpe = 6;
      bType = 9;
      Dim = 3;
      break;
  }
  for(int i=0;i<Ne;i++)
    connect[i] = new int[Npe];
  for(int i=0;i<Nte-Ne;i++)
    bconnect[i] = new int[Nbpe+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-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            break;
          case 2:       //3-point 2D triangle
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            break;
          case 3:       //4-point 2D quadrangle
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            connect[index-Nb][3] = atoi(token[8]);
            break;
          case 4:       //4-point 3D tetrahedra
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            connect[index-Nb][3] = atoi(token[8]);
            break;
          case 8:       //3-point 1D line
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            break;
          case 9:       //6-point 2D triangle
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            connect[index-Nb][3] = atoi(token[8]);
            connect[index-Nb][4] = atoi(token[9]);
            connect[index-Nb][5] = atoi(token[10]);
            break;
          case 10:      //8-point 2D quadrangle
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            connect[index-Nb][3] = atoi(token[8]);
            connect[index-Nb][4] = atoi(token[9]);
            connect[index-Nb][5] = atoi(token[10]);
            connect[index-Nb][6] = atoi(token[11]);
            connect[index-Nb][7] = atoi(token[12]);
            break;
          case 11:      //10-point 3D tetrahedra
            connect[index-Nb][0] = atoi(token[5]);
            connect[index-Nb][1] = atoi(token[6]);
            connect[index-Nb][2] = atoi(token[7]);
            connect[index-Nb][3] = atoi(token[8]);
            connect[index-Nb][4] = atoi(token[9]);
            connect[index-Nb][5] = atoi(token[10]);
            connect[index-Nb][6] = atoi(token[11]);
            connect[index-Nb][7] = atoi(token[12]);
            connect[index-Nb][8] = atoi(token[13]);
            connect[index-Nb][9] = atoi(token[14]);
            break;
          }
        }
        else{
          switch (Type)
          {
          case 1:       //1-point (2-point 1D line)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            break;
          case 2:       //2-point line (3-point 2D triangle)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            bconnect[index][2] = atoi(token[6]);
            break;
          case 3:       //2-point line (4-point 2D quadrangle)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            bconnect[index][2] = atoi(token[6]);
            break;
          case 4:       //3-point triangle (4-point 3D tetrahedra)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            bconnect[index][2] = atoi(token[6]);
            bconnect[index][3] = atoi(token[7]);
            break;
          case 8:       //1-point (3-point 1D line)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            break;
          case 9:       //3-point line (6-point 2D triangle)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            bconnect[index][2] = atoi(token[6]);
            bconnect[index][3] = atoi(token[7]);
            break;
          case 10:      //3-point line (8-point 2D quadrangle)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            bconnect[index][2] = atoi(token[6]);
            bconnect[index][3] = atoi(token[7]);
            break;
          case 11:      //6-point triangle (10-point 3D tetrahedra)
            bconnect[index][0] = atoi(token[3]);
            bconnect[index][1] = atoi(token[5]);
            bconnect[index][2] = atoi(token[6]);
            bconnect[index][3] = atoi(token[7]);
            bconnect[index][4] = atoi(token[8]);
            bconnect[index][5] = atoi(token[9]);
            bconnect[index][6] = atoi(token[10]);
            break;
          }
        }
      }
    }
    index++;
  }
  fin.close();
}  




//----------------------------------------------------------------------------
// Read model.txt file specifying boundary conditions and material properties
//----------------------------------------------------------------------------
void FeModel::readModel()
{
  ifstream gin;
  gin.open("model.txt");

  nDirGrps = 0;
  nHFGrps = 0;
  nCGrps = 0;

  while(!gin.eof())
  {
    char buf[MAX_CHARS_PER_LINE];
    gin.getline(buf,MAX_CHARS_PER_LINE);
    const char* token[MAX_TOKENS_PER_LINE]={};

    int n=0;
    int flag=0;
    token[0] = strtok(buf,DELIMITER);

    while(true){
      if(!token[n]){break;}
      if(strcmp(token[n],"#")==0){break;}
      if(strcmp(token[n],"dt")==0){flag=1;}
      if(strcmp(token[n],"steps")==0){flag=2;}
      if(strcmp(token[n],"E")==0){flag=3;}
      if(strcmp(token[n],"mu")==0){flag=4;}
      if(strcmp(token[n],"alpha")==0){flag=5;}
      if(strcmp(token[n],"T")==0){flag=6;}
      if(strcmp(token[n],"time")==0){flag=7;}
      if(strcmp(token[n],"dirichletGrps")==0){flag=8;}
      if(strcmp(token[n],"dirichletVals")==0){flag=9;}
      if(strcmp(token[n],"heatFluxGrps")==0){flag=10;}
      if(strcmp(token[n],"heatFluxVals")==0){flag=11;}
      if(strcmp(token[n],"convectionGrps")==0){flag=12;}
      if(strcmp(token[n],"convectionVals")==0){flag=13;}
      n++;
      token[n] = strtok(0,DELIMITER);
    }

    switch(flag)
    {
      case 1:
        dt=strtod(token[2],NULL);
        break;
      case 2:
        steps=atoi(token[2]);
        break;
      case 3:
        E=strtod(token[2],NULL);
        break;
      case 4:
        mu=strtod(token[2],NULL);
        break;
      case 5:
        alpha=strtod(token[2],NULL);
        break;
      case 6:
        T=strtod(token[2],NULL);
        break;
      case 7:
        time = false;
        if(strcmp(token[2],"true")==0){time = true;}
        break;
      case 8:
        nDirGrps = n-2;
        dirichletGrps = new int[n-2];
        for(int i=2;i<n;i++){
          dirichletGrps[i-2] = atoi(token[i]);
        }
        break;
      case 9:
        dirichletVals = new double[n-2];
        for(int i=2;i<n;i++){
          dirichletVals[i-2] = strtod(token[i],NULL);
        }
        break;  
      case 10:
        nHFGrps = n-2;
        heatFluxGrps = new int[n-2];
        for(int i=2;i<n;i++){
          heatFluxGrps[i-2] = atoi(token[i]);
        }
        break;
      case 11:
        heatFluxVals = new double[n-2];
        for(int i=2;i<n;i++){
          heatFluxVals[i-2] = strtod(token[i],NULL);
        }
        break;
      case 12:
        nCGrps = n-2;
        convectionGrps = new int[n-2];
        for(int i=2;i<n;i++){
          convectionGrps[i-2] = atoi(token[i]);
        }
        break;
      case 13:
        convectionVals = new double[n-2];
        for(int i=2;i<n;i++){
          convectionVals[i-2] = strtod(token[i],NULL);
        }
        break;
    }
  }

  display();
  gin.close();
}




//----------------------------------------------------------------------------
// Delete Model
//----------------------------------------------------------------------------
void FeModel::deleteModel()
{
  delete [] xpoints;
  delete [] ypoints;
  delete [] zpoints;

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

  for(int i=0;i<Nb;i++){
    delete [] bconnect[i];
  }
  delete [] bconnect;

  delete [] dirichletGrps;
  delete [] dirichletVals;
}



//----------------------------------------------------------------------------
// Display Model Header
//----------------------------------------------------------------------------
void FeModel::display()
{
  cout<<"**************************************************************"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"*          C++ Finite Element Solver for heat equation       *"<<endl;
  cout<<"*                     James Sandham                          *"<<endl;
  cout<<"*                      05 May 2014                           *"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"* This code uses gmsh (.msh) ASCII files for mesh input. As  *"<<endl;
  cout<<"* per the gmsh standard, mesh elements are coded as:         *"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"* Linear Line (1D) --------------------------------- 1       *"<<endl;
  cout<<"* Linear Triangles (2D) ---------------------------- 2       *"<<endl;
  cout<<"* Linear Quadrangles (2D) -------------------------- 3       *"<<endl;
  cout<<"* Linear Tetrahedra (3D) --------------------------- 4       *"<<endl;
  cout<<"* Quadratic Line (1D) ------------------------------ 8       *"<<endl;
  cout<<"* Quadratic Triangles (2D) ------------------------- 9       *"<<endl;
  cout<<"* Quadratic Quadrangles (2D) ----------------------- 10      *"<<endl;
  cout<<"* Quadratic Tetrahedra (3D) ------------------------ 11      *"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"**************************************************************"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"* Np = "<<N<<"                                              *"<<endl;
  cout<<"* Nte = "<<Nte<<"                                            *"<<endl;
  cout<<"* Ne = "<<Ne<<"                                              *"<<endl;
  cout<<"* Nb = "<<Nb<<"                                              *"<<endl;
  cout<<"* Ng = "<<Ng<<"                                              *"<<endl;
  cout<<"* Npe = "<<Npe<<"                                            *"<<endl;
  cout<<"* Nbpe = "<<Nbpe<<"                                          *"<<endl;
  cout<<"* Type = "<<Type<<"                                          *"<<endl;
  cout<<"* bType = "<<bType<<"                                        *"<<endl;
  cout<<"*                                                            *"<<endl;
  cout<<"**************************************************************"<<endl;
}
