// 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 "GmshGlobal.h"

#include "genMatrix.h"
#include "genTensors.h"
#include "genSupports.h"
#include "genSimpleFunction.h"
#include "genSimpleFunctionExtended.h"

#include "ballManager.h"
#include "genFilters.h"

#include "GModel.h"
#include "MElement.h"

class Data
{
public :
  double distance;
  Data(double distance_) : distance(distance_) {}
};

int main (int argc, char* argv[])
{
  GmshInitialize(argc, argv);

/////////////////////////
// STEP 0 - Parameters //
/////////////////////////
//   std::string GModelFileName="rectangle.msh";
//LevelSet 0 1 Function Scalar 0.*x+1.*y+0.*z-1.5 [x,y,z]
//   std::string GModelFileName="rectangle_*3.msh";
//LevelSet 0 1 Function Scalar 0.*x+1.*y+0.*z-5.5 [x,y,z]
  std::string GModelFileName="rve_cut.msh";
//LevelSet 1 2 Function Scalar sqrt(((0-x)/0.4)^2+((0-y)/0.4)^2)-1 [x,y,z]
  int tag=1, dim=2;
  std::vector<std::string> expr(1),var(3);
  expr[0]="sqrt(((0-x)/0.4)^2+((0-y)/0.4)^2)-1";
  var[0]="x"; var[1]="y"; var[2]="z";
  genSimpleFunction<genTensor0<double> >* fscalar = new genSimpleFunctionAnalytical<genTensor0<double> >(expr,var);
  std::map<int,std::vector<genScalar<genTensor0<double> > > > distances; // vecteur des distances des noeuds de l element, respectant l'ordre interne de l'element


//////////////////////////////////////
// STEP 1 - Build Group of Elements //
//////////////////////////////////////
  GModel* Model = new GModel();
  int readingDone = Model->readMSH(GModelFileName.c_str());
  assert(readingDone);
  BallManager<Data,int>* ballManager = BallManager<Data,int>::getInstance();

//--------------------------------------------------------------------------
  printf("---> Start bulding ballManager\n");
  genGroupOfElements* allElementsWithoutChildren = new genGroupOfElements();
  genFilterElementParent filterParent;

  // Elements parents 3D
  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);

  // Elements parents 2D
  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);

  // Elements parents 1D
  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");


/////////////////////////////////
// STEP 2 - Build Distance map //
/////////////////////////////////

  // génère les données : distances, (normales et courbure).
  printf("   ---> the corresponding genGroupOfElements already exists\n");
  genGroupOfElements ge(dim, tag);
  genFilterElementTrivial filter;
  genGroupOfElements gp(ge.pbegin(), ge.pend(), filter);
  genGroupOfElements* 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;
  }


////////////////////////////
// STEP 3 - Print results //
////////////////////////////
  printf("\n#### BElements ####\n");
  BallManager<Data,int>::BElementContainer::const_iterator itBe;
  for (itBe = ballManager->begin(); itBe != ballManager->end(); ++itBe)
  {
    BElement<int,Data>* Be = (*itBe);
    std::cout << (*Be);
  }

  printf("\n#### BVertices ####\n");
  BallManager<Data,int>::BVertexContainer::const_iterator itBv;
  for (itBv = ballManager->vbegin(); itBv != ballManager->vend(); ++itBv)
  {
    BVertex<Data,int>* Bv = (*itBv);
    std::cout << (*Bv);
    if (Bv->flag) std::cout << " distance " << Bv->getData()->distance << std::endl;
  }

  printf("\n#### Results ####\n");
  std::map<int,std::vector<genScalar<genTensor0<double> > > >::iterator itd;
  for (itd=distances.begin(); itd!=distances.end();++itd)
  {
    printf("Element[%d] = [",itd->first);
    for (int i=0; i<itd->second.size(); ++i)
    {
      printf("%lf",itd->second[i]()());
      if (i+1<itd->second.size()) printf(", ");
    }
    printf("]\n");
  }

  // stop gmsh
  GmshFinalize();
}