// gen_LevelSet - An automatic generator of meshes from level-sets
// Copyright (C) 2012-2026 Eric Bechet, Frederic Duboeuf
//
// See the LICENSE file for license information.
// Please report all bugs and problems to <bechet@cadxfem.org> or <duboeuf@outlook.com>.


#include "genDataLS.h"
#include "genFilters.h"
#include "ballManager.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cassert>
#include "GModel.h"
#include "genSimpleFunctionExtended.h"

class MElement;

void genDataLS::readToken(std::istream &is, std::string &token)
{
  if (token=="LevelSet")
  {
    genLevelSetSupport *LS = new genLevelSetSupport;
    is >> *LS;
    LS->group();
    LSs.push_back(LS);
  }
  else
  {
    genDataIO::readToken(is, token);
  }
}

void genDataLS::dumpContent(std::ostream &os)
{
  genDataIO::dumpContent(os);
  for (int i=0;i<LSs.size();++i) os << *LSs[i] << std::endl;
}


class Data
{
public :
  double distance;
  Data(double distance_) {distance = (fabs(distance_)<TOL) ? 0. : distance_;}
};

genGroupOfElements* genLevelSetSupport::group()
{
  if (g==NULL)
  {
    BallManager<Data,int>* ballManager = BallManager<Data,int>::getInstance();
    if (!ballManager->isItAlreadyBuilt())
    {
      printf("---> Start bulding ballManager\n");
      genGroupOfElements* allElementsWithoutChildren = new genGroupOfElements();
      genFilterElementParent filterParent;

      GModel::riter itGR;
      for (itGR = GModel::current()->firstRegion(); itGR != GModel::current()->lastRegion(); ++itGR)
      {
        GRegion* gr = *itGR;
        genFilterElementType filterType((int)TYPE_TET);
        genFilterElementOperatorAnd filter(&filterParent,&filterType);
        allElementsWithoutChildren->addElementary(gr, filter);
      }
      int nbtetra = allElementsWithoutChildren->size();
      printf("   ---> genGroupOfElements : %d tetrahedron elements\n", nbtetra);

      GModel::fiter itGF;
      for (itGF = GModel::current()->firstFace(); itGF != GModel::current()->lastFace(); ++itGF)
      {
        GFace* gf = *itGF;
        genFilterElementType filterType((int)TYPE_TRI);
        genFilterElementOperatorAnd filter(&filterParent,&filterType);
        allElementsWithoutChildren->addElementary(gf, filter);
      }
      int nbtri = allElementsWithoutChildren->size()-nbtetra;
      printf("   ---> genGroupOfElements : %d triangle elements\n", nbtri);

      GModel::eiter itGE;
      for (itGE = GModel::current()->firstEdge(); itGE != GModel::current()->lastEdge(); ++itGE)
      {
        GEdge* ge = *itGE;
        genFilterElementType filterType((int)TYPE_LIN);
        genFilterElementOperatorAnd filter(&filterParent,&filterType);
        allElementsWithoutChildren->addElementary(ge, filter);
      }
      int nblin = allElementsWithoutChildren->size()-(nbtetra+nbtri);
      printf("   ---> genGroupOfElements : %d line elements\n", nblin);

      printf("   ---> genGroupOfElements : %d vertices, %d elements\n", allElementsWithoutChildren->vsize(), allElementsWithoutChildren->size());
      ballManager->setGroup(allElementsWithoutChildren);
      printf("   ---> new ballManager   : %d vertices, %d elements\n", ballManager->vsize(), ballManager->size());
      printf("---> End bulding ballManager\n");
    }
    else
    {
      ballManager->clearDatas();
      printf("---> Use existing ballManager\n");
    }

    printf("---> Start bulding genLevelSetSupport : dim=%d, tag=%d\n", dim, tag);

    if (kind=="Mesh")
    {// construit les données : distances, (normales et courbure).
      // sélectionne les éléments du support en leur associant une distance
    }
    else if (kind=="Function")
    {// génère les données : distances, (normales et courbure).
      printf("   ---> use Function to compute distances\n");

      distances = new std::map<int,std::vector<genScalar<genTensor0<double> > > >;

      // cherche s'il existe déjà le groupe avec le physical tag de dimension dim
      std::map<int, std::vector<GEntity*> > groups[4];
      GModel::current()->getPhysicalGroups(groups);
      std::map<int, std::vector<GEntity*> >::iterator ite = groups[dim].find(tag);

      if (ite==groups[dim].end())
      { // si le tag n'existe pas déjà
        // alors on génère le support et les distances
        printf("   ---> the corresponding genGroupOfElements does not already exist\n");
        meshContainsLS = 0;
        g = new genGroupOfElements;

        // On remplie les données relatives aux distances
        BallManager<Data,int>::BVertexContainer::const_iterator itBv;
        for (itBv = ballManager->vbegin(); itBv != ballManager->vend(); ++itBv)
        {
          BallManager<Data,int>::BVertexType* Bv = *itBv;
          MVertex* v = Bv->v;
          assert(!Bv->getData());
          double distance = (*fscalar)( v->x(), v->y(), v->z() );
          Bv->setData( new Data(distance) );
        }

        // Pour ajouter un élément dans le support :
        // soit l'interface intersecte l'élément
        // soit l'un des sommets a déjà été identifié comme support
        // sinon ne rien faire

        // On sélectionne les sommets du support
        std::vector<BallManager<Data,int>::BVertexType*> selectionBv;
        BallManager<Data,int>::BElementContainer::const_iterator itBe;
        for (itBe = ballManager->begin(); itBe != ballManager->end(); ++itBe)
        {
          BallManager<Data,int>::BElementType* Be = *itBe;

          double d[Be->getNumBvertices()];

          for (int j=0; j<Be->getNumBvertices(); ++j)
            d[j] = Be->getBVertice(j)->getData()->distance;

          bool addToSelection = false;
          // l'interface intersecte elle l'élément ?
          switch (Be->e->getType())
          {
            case TYPE_LIN:
              if (d[0]*d[1]<=0) addToSelection = true;
              break;
            case TYPE_TRI:
              if (d[0]*d[1]<=0) addToSelection = true;
              else if (d[0]*d[2]<=0) addToSelection = true;
              else if (d[1]*d[2]<=0) addToSelection = true;
              break;
            case TYPE_TET:
              if (d[0]*d[1]<=0) addToSelection = true;
              else if (d[0]*d[2]<=0) addToSelection = true;
              else if (d[0]*d[3]<=0) addToSelection = true;
              else if (d[1]*d[2]<=0) addToSelection = true;
              else if (d[1]*d[3]<=0) addToSelection = true;
              else if (d[2]*d[3]<=0) addToSelection = true;
              break;
          }

          if (addToSelection)
          {
            for (int j=0; j<Be->getNumBvertices(); ++j)
            {
              if (!Be->getBVertice(j)->flag)
              {
                Be->getBVertice(j)->flag = 1;
                selectionBv.push_back(Be->getBVertice(j));
              }
            }
          }
        }
        printf("   ---> %d vertices selected\n", (int)selectionBv.size());

        // On sélectionne les éléments du support et de son plus proche voisinage
        std::vector<BallManager<Data,int>::BElementType*> selectionBe;
        for (int i=0; i<selectionBv.size(); ++i)
        {
          BallManager<Data,int>::BVertexType* Bv = selectionBv[i];
          for (int j=0; j<Bv->getNumBelements(); ++j)
          {
            BallManager<Data,int>::BElementType* Be = Bv->getBElement(j);
            if (!Be->flag)
            {
              Be->flag = 1;
              selectionBe.push_back(Bv->getBElement(j));
            }
          }
        }
        printf("   ---> %d elements selected\n", (int)selectionBe.size());

        // On construit le support et les distances associées
        for (int i=0; i<selectionBe.size(); ++i)
        {
          BallManager<Data,int>::BElementType* Be = selectionBe[i];
          g->insert(Be->e);
          std::vector<genScalar<genTensor0<double> > > d(Be->getNumBvertices());
          for (int j=0; j<Be->getNumBvertices(); ++j)
          {
            BallManager<Data,int>::BVertexType* Bv = Be->getBVertice(j);
            d[j] = genScalar<genTensor0<double> >(Bv->getData()->distance);
          }
          (*distances)[Be->getNum()] = d;
        }
      }

      else
      { // si le tag existe déjà
        // alors on charge simplement le support et l'on génère les distances
        printf("   ---> the corresponding genGroupOfElements already exists\n");
        meshContainsLS = 1;
        genGroupOfElements ge(dim, tag);
        genFilterElementTrivial filter;
        genGroupOfElements gp(ge.pbegin(), ge.pend(), filter);
        g = new genGroupOfElements;
        // utiliser ce groupe comme éléments de départ intersectant le support
        // récupérer les éléments du plus proche voisinage
        // puis calculer les distances pour chacun des sommets
        // pour enfin construire les elements et distances

        // On sélectionne les éléments du support et de son plus proche voisinage
        std::vector<BallManager<Data,int>::BElementType*> selectionBe;
        std::set<MVertex*>::const_iterator itv;
        for (itv = ge.vbegin(); itv != ge.vend(); ++itv)
        {
          MVertex* v = *itv;
          BallManager<Data,int>::BVertexContainer::iterator itBv;
          itBv = ballManager->vfind(v);
          assert(itBv!=ballManager->vend());
          BallManager<Data,int>::BVertexType* Bv = *itBv;

          for (int j=0; j<Bv->getNumBelements(); ++j)
          {
            BallManager<Data,int>::BElementType* Be = Bv->getBElement(j);
            if (!Be->flag)
            {
              Be->flag = 1;
              selectionBe.push_back(Bv->getBElement(j));
            }
          }
        }
        printf("   ---> %d vertices selected\n", (int)ge.vsize());
        printf("   ---> %d elements selected\n", (int)selectionBe.size());

        // On construit le support et les distances associées
        for (int i=0; i<selectionBe.size(); ++i)
        {
          BallManager<Data,int>::BElementType* Be = selectionBe[i];
          g->insert(Be->e);
          std::vector<genScalar<genTensor0<double> > > d(Be->getNumBvertices());
          for (int j=0; j<Be->getNumBvertices(); ++j)
          {
            BallManager<Data,int>::BVertexType* Bv = Be->getBVertice(j);
            if (!Bv->flag)
            {
              MVertex* v = Bv->v;
              assert(!Bv->getData());
              double distance = (*fscalar)( v->x(), v->y(), v->z() );
              Bv->setData( new Data(distance) );
              Bv->flag = 1;
            }
            d[j] = genScalar<genTensor0<double> >(Bv->getData()->distance);
          }
          (*distances)[Be->getNum()]=d;
        }
      }
    }
    else if (kind=="File")
    {// lit les données : distances, (normales et courbure).
      // importe les éléments du support avec leurs distances associées
    }
    printf("---> End bulding genLevelSetSupport : %d vertices, %d elements\n", g->vsize(), g->size());
    printf("---> End bulding distances : %d elements\n", (int)distances->size());
    assert(g->size()!=0);
  }
  return g;
}

std::istream & operator>>(std::istream &is, genLevelSetSupport &obj)
{
  int physical;
  int dim;
  readScalar(is,physical);
  readScalar(is,dim);
  obj.tag = physical;
  obj.dim = dim;
  NextTokenWithoutComments(is,obj.kind);
  if (obj.kind=="Mesh")
  {
    // on doit construire le support (elt parent) et non uniquement les elements de l'interface.
    // obj.distances ...
  }
  else if (obj.kind=="Function")
  {
    std::string iswhat;
    NextTokenWithoutComments(is,iswhat); // Scalar,Vector
    assert (iswhat=="Scalar");
    {
      std::vector<std::string> expr(1),var;
      readString(is,expr[0]);
      readVector(is,var);
      obj.fscalar = new genSimpleFunctionAnalytical<genTensor0<double> >(expr,var);
    }
  }
  else if (obj.kind=="File")
  {
    NextTokenWithoutComments(is,obj.LSfileName);
  }
  return is;
}

std::ostream & operator<<(std::ostream &os, const genLevelSetSupport &obj)
{
  os << "LevelSet " << obj.tag << " " << obj.dim << " " << obj.kind;
  if (obj.kind=="File") os << obj.LSfileName << " ";
  if (obj.kind=="Function") os << obj.LSfileName << " " << "Scalar ";
  return os;
}

void writeLSs(std::vector<genLevelSetSupport*> &LSs, const std::string &baseName)
{
  char lsfileName[100];
  sprintf(lsfileName, "%s.ls", baseName.c_str());
  std::ofstream ofs;
  ofs.open(lsfileName);

  ofs << "0" << "\n"; // method
  ofs << (int)LSs.size() << "\n"; // nb_LS

  for (int i=0;i<LSs.size();++i)
  {
    genLevelSetSupport* LS = LSs[i];
    writeLS(ofs, LS, i);
  }
  ofs << "0" << "\n";
  ofs.close();
}

void writeLS(std::ostream &os, genLevelSetSupport* LS, int LSNum)
{
  int groupSize = LS->group()->size();
  assert(groupSize==LS->distances->size());

  os << LSNum << " " << groupSize << "\n";

  std::set<MElement*>::const_iterator ite=LS->begin();
  std::map<int,std::vector<genScalar<genTensor0<double> > > >::const_iterator itd;
  for (ite = LS->begin(); ite != LS->end(); ++ite)
  {
    MElement* e = *ite;

    int eNum = e->getNum();
    itd = LS->distances->find(eNum);
    assert(itd!=LS->distances->end());

    os << eNum << " ";

    std::vector<MVertex*> v;
    e->getVertices(v);
    for (int k=0; k<(int)v.size(); ++k)
    {
      os << v[k]->getNum() << " " << (double)itd->second[k]() << " ";
    }
    os << "\n";
  }
  os << "\n";
}