// GenFem - A high-level finite element library
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//



#include "genSimpleFunctionExtended.h"
#include "GModel.h"
#include "genDataIO.h"

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cassert>
#include <boost/algorithm/string.hpp>

//const char genBoundaryConditionNG::physNames[][256]={"Displacement","Force","Potential","Charge","Temperature","ThermalFlux"};

genBoundaryConditionNG::~genBoundaryConditionNG()
{
  if (own)
  {
    if (fscalar) delete fscalar;
    if (fvector) delete fvector;
    if (ftensor) delete ftensor;
    if (fcscalar) delete fcscalar;
    if (fcvector) delete fcvector;
    if (fctensor) delete fctensor;
    if (periodic) delete periodic;
  }
}


genDataIONG::~genDataIONG()
{
  for (int i=0;i<Domains.size();++i) delete Domains[i];
  for (int i=0;i<BCs.size();++i) delete BCs[i];
  if (Model) delete Model;
}

void genDataIONG::readInputFile(const std::string &fileName)
{
  std::ifstream ifs(fileName.c_str());
  if (!ifs)
  {
    printf("Unable to open file '%s'\n", fileName.c_str());
    exit(0);
  }
  while ( !ifs.eof() )
  {
    std::string token;
    NextTokenWithoutComments(ifs, token);
    readToken(ifs, token);
  }
  ifs.close();
}

void genDataIONG::readToken(std::istream &is, std::string &token)
{
  if ( token=="Dim")
  {
    readScalar(is,dim);
  }
  else if (token=="MeshFile")
  {
    readString(is,GModelBaseName);
    std::string::size_type idx = GModelBaseName.find(".msh");
    GModelBaseName = GModelBaseName.substr(0, idx).c_str();

    std::string GModelFileName;
    GModelFileName = GModelBaseName + ".msh";

    printf("mesh name : %s\n", GModelFileName.c_str());
    Model = new GModel();
    int readingDone;
    readingDone = Model->readMSH(GModelFileName.c_str());
    if (!readingDone) exit(0);

    if (dim!=0) dim = dim;
    else if (Model->getNumRegions()) dim = 3;
    else if (Model->getNumFaces())   dim = 2;
    else if (Model->getNumEdges())   dim = 1;
  }
  else if (token=="BoundaryCondition")
  {
    genBoundaryConditionNG* BC = new genBoundaryConditionNG;
    is >> *BC;
    BCs.push_back ( BC );
  }
  else if (token=="Domain" )
  {
    genDomain* Sup=new genDomain;
    is >> *Sup;
    Domains.push_back(Sup);
  }
  else if ( token=="Solver")
  {
    readScalar(is,solvertype);
  }
  else if (token.find_first_of("/#.;,\\+-*[]")==std::string::npos)
  {
    std::string tenstype;
    NextTokenWithoutComments(is,tenstype);
    if (tenstype=="Integer")
    {
      readScalar(is,User.Integers[token]);
      User.AllData[token]=genDomain::Integer;
    }
    else if (tenstype=="Scalar")
    {
      readScalar(is,User.Scalars[token]);
      User.AllData[token]=genDomain::Scalar;
    }
    else if (tenstype=="CScalar")
    {
      readScalar(is,User.CScalars[token]);
      User.AllData[token]=genDomain::CScalar;
    }
    else if (tenstype=="Vector")
    {
      readVector(is,User.Vectors[token]);
      User.AllData[token]=genDomain::Vector;
    }
    else if (tenstype=="CVector")
    {
      readVector(is,User.CVectors[token]);
      User.AllData[token]=genDomain::CVector;
    }
    else if (tenstype=="Matrix")
    {
      readMatrix(is,User.Tensors[token]);
      User.AllData[token]=genDomain::Tensor;
    }
    else if (tenstype=="CMatrix")
    {
      readMatrix(is,User.CTensors[token]);
      User.AllData[token]=genDomain::CTensor;
    }

  }
}



void genDataIONG::dumpContent(std::ostream &os)
{
  if (Model)
  {
    os << "MeshFile " << GModelBaseName << std::endl;
  }
  for (int i=0;i<Domains.size();++i) os << *Domains[i] << std::endl;
  for (int i=0;i<BCs.size();++i) os << *BCs[i] << std::endl;
  for (std::map<std::string,genDomain::tenstype>::const_iterator it=User.AllData.begin();it!=User.AllData.end();++it)
  {
    if (it->second==genDomain::Integer)
      os << it->first << " Integer " << User.Integers.find(it->first)->second << std::endl;
    else if (it->second==genDomain::Scalar)
      os << it->first << " Scalar " << User.Scalars.find(it->first)->second << std::endl;
    else if (it->second==genDomain::CScalar)
      os << it->first << " CScalar " << User.CScalars.find(it->first)->second << std::endl;
    else if (it->second==genDomain::Vector)
      os << it->first << " Vector " << "XXXXXXXXX"/*User.Vectors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::CVector)
      os << it->first << " CVector " << "XXXXXXXXX"/*User.CVectors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::Tensor)
      os << it->first << " Matrix " << "XXXXXXXXX"/*User.Tensors.find(it->first)->second */<< std::endl;
    else if (it->second==genDomain::CTensor)
      os << it->first << " CMatrix " << "XXXXXXXXX"/*User.CTensors.find(it->first)->second */<< std::endl;
  }
  os << "SolverType " << solvertype << std::endl;
}


std::istream & operator>>(std::istream &is, genBoundaryConditionNG &obj)
{
  std::string topo;
  obj.dim=genBoundaryConditionNG::UNDEF;
  obj.kind="";
  do
  {
    NextTokenWithoutComments(is,topo);
    if (topo=="Node") obj.dim=genBoundaryConditionNG::ON_VERTEX;
    else if (topo=="Edge") obj.dim=genBoundaryConditionNG::ON_EDGE;
    else if (topo=="Face") obj.dim=genBoundaryConditionNG::ON_FACE;
    else if (topo=="Volume") obj.dim=genBoundaryConditionNG::ON_VOLUME;
    else obj.kind=topo;
  } while ((obj.dim==genBoundaryConditionNG::UNDEF)&&!is.eof());
  readScalar(is,obj.tag);
  std::string field;
//  obj.What =genBoundaryConditionNG::UserDefined;
  NextTokenWithoutComments(is,obj.What); // DISPLACEMENT , POTENTIAL, etc...
//  for (int i=0;i<(int)genBoundaryConditionNG::Periodic;++i)
//    if (obj.user == genBoundaryConditionNG::physNames[i] ) {obj.What=(genBoundaryConditionNG::physical)i;break;}
  std::string nnn;
  std::string iswhat;
  bool ok=false;
  do
  {
    NextTokenWithoutComments(is,nnn);
    if (nnn=="Constant")
    {
      NextTokenWithoutComments(is,iswhat); // Scalar,Vector
      if (iswhat=="Scalar")
      {
        obj.isWhat=genBoundaryConditionNG::Scalar;
        double val;
        readScalar(is,val);
        obj.fscalar = new genSimpleFunctionConstant<genTensor<0,double> > ( val );
      }
      if (iswhat=="CScalar")
      {
        obj.isWhat=genBoundaryConditionNG::CScalar;
        std::complex<double> val;
        readScalar(is,val);
        obj.fcscalar = new genSimpleFunctionConstant<genTensor<0,std::complex<double> > > ( val );
      }
      else if (iswhat=="Vector")
      {
        obj.isWhat=genBoundaryConditionNG::Vector;
        std::vector<double> val;
        readVector(is,val);
        const double tab[3]={val[0], val[1], val[2]};
        obj.fvector = new genSimpleFunctionConstant<genTensor<1,double> > ( genTensor<1,double> ( &tab[0] ) );
      }
      else if (iswhat=="CVector")
      {
        obj.isWhat=genBoundaryConditionNG::CVector;
        std::vector<std::complex<double> > val;
        readVector(is,val);
        
        obj.fcvector = new genSimpleFunctionConstant<genTensor<1,std::complex<double> > > ( genTensor<1,std::complex<double> > ( &val[0]) );
      }
      else if (iswhat=="Matrix")
      {
        obj.isWhat=genBoundaryConditionNG::Tensor;
        genMatrix<double> val;
        readMatrix(is,val);
        genTensor<2,double> M(val.getDataPtr());
        obj.ftensor = new genSimpleFunctionConstant<genTensor<2,double> > ( M );
      }
      else if (iswhat=="CMatrix")
      {
        obj.isWhat=genBoundaryConditionNG::CTensor;
        genMatrix<std::complex<double> > val;
        readMatrix(is,val);
        genTensor<2,std::complex<double> > M(val.getDataPtr());
        obj.fctensor = new genSimpleFunctionConstant<genTensor<2,std::complex<double> > > ( M );
      }
      ok=false;
    }
    else
    if (nnn=="Function")
    {
      NextTokenWithoutComments(is,iswhat); // Scalar,Vector
      if (iswhat=="Scalar")
      {
          obj.isWhat=genBoundaryConditionNG::Scalar;
          std::vector<std::string> expr(1),var;
          readString(is,expr[0]);
          readVector(is,var);
          obj.fscalar = new genSimpleFunctionAnalytical<genTensor<0,double> >(expr,var);
      }
      if (iswhat=="CScalar")
      {
          obj.isWhat=genBoundaryConditionNG::CScalar;
          std::vector<std::string> expr(1),var;
          readString(is,expr[0]);
          readVector(is,var);
//          obj.fcscalar = new genSimpleFunctionAnalytical<genTensor<0,std::complex<double> > >(expr,var);
      }
      else if (iswhat=="Vector")
      {
        obj.isWhat=genBoundaryConditionNG::Vector;
        std::vector<std::string> expr,var;
        readVector(is,expr);
        readVector(is,var);
        obj.fvector = new genSimpleFunctionAnalytical<genTensor<1,double> >(expr,var);
      }
      else if (iswhat=="CVector")
      {
        obj.isWhat=genBoundaryConditionNG::CVector;
        std::vector<std::string> expr,var;
        readVector(is,expr);
        readVector(is,var);
//        obj.fcvector = new genSimpleFunctionAnalytical<genTensor<1,std::complex<double> > >(expr,var);
      }
      ok=false;
    }
    else
    if (nnn=="Periodic")
    {
      obj.kind="Periodic";
      obj.periodic=new genBoundaryConditionNG;
      is>>*(obj.periodic);
      std::vector<std::string> expr,var;
      readVector(is,expr);
      readVector(is,var);
      obj.fvector = new genSimpleFunctionAnalytical<genTensor<1,double> >(expr,var);
      ok=false;
    }
    else
    {
      obj.comp1=atoi(nnn.c_str());
      ok=true;
    }
  } while (ok);
  return is;
}

std::ostream & operator<<(std::ostream &os, const genBoundaryConditionNG &obj)
{
  os << "BoundaryCondition ";
  if (obj.dim==genBoundaryConditionNG::ON_VERTEX) os << "Vertex ";
  else if (obj.dim==genBoundaryConditionNG::ON_EDGE) os << "Edge ";
  else if (obj.dim==genBoundaryConditionNG::ON_FACE) os << "Face ";
  else if (obj.dim==genBoundaryConditionNG::ON_VOLUME) os << "Volume ";
  os << obj.tag << " " << obj.isWhat ;
  if (obj.periodic!=NULL) os << " Periodic ";

  return os;
}