// read_ucd.cc   read an AVS UCD file, typically  .inp
//               static data in read_ucd.h

#include <iostream>
#include <fstream>
#include <vector>
#include <cctype>
#include "read_ucd.h"

using namespace std;
static bool debug = false;

void skip_comment_lines (std::istream &in,
			 const char    comment_start)
{
  char c, cb;
  // loop over the following comment
  // lines
  while ((c=in.get()) == comment_start)
    // loop over the characters after
    // the comment starter
    while ((cb=in.get()) != '\n');
  
  
  // put back first character of
  // first non-comment line
  in.putback (c);

} // end skip_comment_lines

void read_ucd (std::istream &in)
{
  // skip comments at start of file
  skip_comment_lines (in, '#');

  if(debug) cout << "read_ucd called" << endl;

  // UCD1 lengths
  in >> n_vertices
     >> n_cells
     >> n_ndata         // number of data vectors
     >> n_cdata         // cell data
     >> n_mdata;        // model data
  in.ignore(256,'\n');

  if(debug) cout << "n_vertices = " << n_vertices << endl;
  if(debug) cout << "n_cells =    " << n_cells << endl;
  if(debug) cout << "n_ndata =    " << n_ndata << endl;
  if(debug) cout << "n_cdata =    " << n_cdata << endl;
  if(debug) cout << "n_mdata =    " << n_mdata << endl;

  // set up array of vertices
  vertices.reserve(n_vertices+1);
  if(debug) cout << "vertices.capacity= " << vertices.capacity() << endl;  

  // UCD2 nodes
  for(int vertex=0; vertex<n_vertices; ++vertex) 
  {
    int vertex_number;
    double x, y, z; 

    // read vertex
    in >> vertex_number >> x >> y >> z;
    in.ignore(256,'\n');
    if(vertex_number<1 || vertex_number>n_vertices)
    { 
      cout << "bad vertex " << vertex_number << ", x=" << x 
	   << ", y=" << y << ", z=" << z << endl;
    }
    else
    {
      if(debug) cout << vertex_number << ", x=" << x << ", y=" << y 
                     << ", z=" << z << endl; 
      // store vertex
      vertices[vertex_number].x = x;
      vertices[vertex_number].y = y;
      vertices[vertex_number].z = z;
    }
  } // end UCD2

  // set up array of cells
  cells.reserve(n_cells+1);
  if(debug) cout << "cells.capacity= " << cells.capacity() << endl;  

  // UCD3 cells
  for (int i=0; i<n_cells; ++i) 
  {
    std::string cell_type;
    int cell_num, node_num, material_id, celli;
      
    in >> cell_num        // cell number
       >> material_id;
    in >> cell_type;
    if(debug) cout << "cell num= " << cell_num << ", matl=" << material_id 
                   << ", type=" << cell_type << endl; 

    if(cell_num<1 || cell_num>n_cells)
    { 
      cout << "bad cell " << cell_num << endl;
      continue;
    } 
    celli = cell_num;

    if (cell_type == "pt")         // 1 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 1;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
    }
    else if (cell_type == "line")  // 2 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 2;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
    }
    else if (cell_type == "tri")   // 3 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 3;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
      in >> node_num;
      cells[celli].nodes[2] = node_num;
    }
    else if (cell_type == "quad")  // 4 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 4;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
      in >> node_num;
      cells[celli].nodes[2] = node_num;
      in >> node_num;
      cells[celli].nodes[3] = node_num;
    }
    else if (cell_type == "tet")   // 4 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 5;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
      in >> node_num;
      cells[celli].nodes[2] = node_num;
      in >> node_num;
      cells[celli].nodes[3] = node_num;
    }
    else if (cell_type == "pyr")   // 5 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 6;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
      in >> node_num;
      cells[celli].nodes[2] = node_num;
      in >> node_num;
      cells[celli].nodes[3] = node_num;
      in >> node_num;
      cells[celli].nodes[4] = node_num;
    }
    else if (cell_type == "prism") // 6 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 7;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
      in >> node_num;
      cells[celli].nodes[2] = node_num;
      in >> node_num;
      cells[celli].nodes[3] = node_num;
      in >> node_num;
      cells[celli].nodes[4] = node_num;
      in >> node_num;
      cells[celli].nodes[5] = node_num;
    }
    else if (cell_type == "hex" )  // 8 vertex
    {
      cells[celli].matl = material_id;
      cells[celli].cell_type = 8;
      in >> node_num;
      cells[celli].nodes[0] = node_num;
      in >> node_num;
      cells[celli].nodes[1] = node_num;
      in >> node_num;
      cells[celli].nodes[2] = node_num;
      in >> node_num;
      cells[celli].nodes[3] = node_num;
      in >> node_num;
      cells[celli].nodes[4] = node_num;
      in >> node_num;
      cells[celli].nodes[5] = node_num;
      in >> node_num;
      cells[celli].nodes[6] = node_num;
      in >> node_num;
      cells[celli].nodes[7] = node_num;
    }
    else
    {
      cout << "bad cell, unknown type is " << cell_type << endl;
    }
    in.ignore(256,'\n');
  } // end UCD3

  // UCD4 node data
  if(n_ndata>0)
  {
    char snodes1[255];
    char snodes2[255];
    char snodes3[255];
    int k;
    node_data.reserve(n_vertices+1);
    if(debug) cout << "node_data.capacity= " << node_data.capacity() << endl;  
    in >> nodecomp >> noden1;
    if(debug) cout << "node components " << nodecomp << " " << noden1;
    if(nodecomp>1) {in >> noden2; if(debug) cout << " " << noden2;}
    if(nodecomp>2) {in >> noden3; if(debug) cout << " " << noden3;}
    in.ignore(256,'\n'); // max 3 implemented
    if(debug) cout << endl;
    in.getline(snodes1,255,'\n');
    nodes1=string(snodes1);
    if(debug) cout << nodes1 << endl;
    if(nodecomp>1) {in.getline(snodes2,255,'\n'); nodes2=string(snodes2);
                    if(debug) cout << nodes2 << endl;}
    if(nodecomp>2) {in.getline(snodes3,255,'\n'); nodes3=string(snodes3);
                    if(debug) cout << nodes3 << endl;}
    if(nodecomp>3) {cout << "too many node components, max=3" << endl;}
    for(int i=0; i<n_vertices; i++)
    {
      in >> k; // node number
      if(debug) cout << "node num=" << k << endl;
      for(int j=0; j<n_ndata; j++)
      {
        in >> node_data[k].val[j];
        if(debug) cout << "j=" << j << " val=" << node_data[k].val[j] << endl;
      }
      in.ignore(256,'\n');
    }   
  } // end UCD4

  // UCD5 cell data
  if(n_cdata>0)
  {
    char scells1[255];
    char scells2[255];
    char scells3[255];
    int k;
    cell_data.reserve(n_cells+1);
    if(debug) cout << "cell_data.capacity= " << cell_data.capacity() << endl;  
    in >> cellcomp >> celln1;
    if(debug) cout << "cell components " << cellcomp << " " << celln1;
    if(cellcomp>1) {in >> celln2; cout << " " << celln2;}
    if(cellcomp>2) {in >> celln3; cout << " " << celln3;}
    in.ignore(256,'\n'); // max 3 implemented
    if(debug) cout << endl;
    in.getline(scells1,255,'\n');
    cells1=string(scells1);
    if(debug) cout << cells1 << endl;
    if(cellcomp>1) {in.getline(scells2,255,'\n'); cells2=string(scells2);
                    if(debug) cout << cells2 << endl;}
    if(cellcomp>2) {in.getline(scells3,255,'\n'); cells3=string(scells3);
                    if(debug) cout << cells3 << endl;}
    if(cellcomp>3) {cout << "too many cell components, max=3" << endl;}
    for(int i=0; i<n_cells; i++)
    {
      in >> k; // cell number
      if(debug) cout << "cell num=" << k << endl;
      for(int j=0; j<n_cdata; j++)
      {
        in >> cell_data[k].val[j];
        if(debug) cout << "j=" << j << " val=" << cell_data[k].val[j] << endl;
      }
      in.ignore(256,'\n');
    }   
  } // end UCD5

  // UCD5 model data
  if(n_mdata>0)
  {
    char smodels1[255];
    char smodels2[255];
    char smodels3[255];
    int k, numcomp;
    model_data.reserve(n_mdata+1);
    if(debug) cout << "model_data.capacity= " << model_data.capacity() << endl;  
    in >> modelcomp >> modeln1;
    numcomp = modeln1;
    if(debug) cout << "model components " << modelcomp << " " << modeln1;
    if(modelcomp>1) {in >> modeln2; numcomp += modeln2; 
                     if(debug) cout << " " << modeln2;}
    if(modelcomp>2) {in >> modeln3; numcomp += modeln3;
                     if(debug) cout << " " << modeln3;}
    in.ignore(256,'\n'); // max 3 implemented
    if(debug) cout << endl;
    in.getline(smodels1,255,'\n');
    models1=string(smodels1);
    if(debug) cout << models1 << endl;
    if(modelcomp>1) {in.getline(smodels2,255,'\n'); models2=string(smodels2);
                    if(debug) cout << models2 << endl;}
    if(modelcomp>2) {in.getline(smodels3,255,'\n'); models3=string(smodels3);
                    if(debug) cout << models3 << endl;}
    if(modelcomp>3) {cout << "too many model components, max=3" << endl;}
    for(int i=0; i<n_mdata; i++)
    {
      in >> k; // model number
      if(debug) cout << "model num=" << k << endl;
      for(int j=0; j<numcomp; j++)
      {
        in >> model_data[k].val[j];
        if(debug) cout << "j=" << j << " val=" << model_data[k].val[j] << endl;
      }
      in.ignore(256,'\n');
    }   
  } // end UCD6

} // end read_ucd