// 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 <cstring>
#include <cstdio>
#include <fstream>

#include "genLevelSet.h"

#include "assert.h"


GModel* genLevelSet::process(GModel* model_, std::vector<LevelSetDomain*> &LSDomains, const std::string &fileName)
{
  model = model_;

  std::vector<LevelSetDomain*> LSforCutting;
  LSforCutting = selectLSforCutting(LSDomains);

  if (LSforCutting.size())
  {
    printf("--> Cut the mesh with %d LSs \n", (int)LSforCutting.size());

    // mesh cutting
    cutModel = cutmesh(model, LSforCutting);

    if (fileName!="") writeCutModel(fileName);
  }
  return cutModel;
}


std::vector<LevelSetDomain*> genLevelSet::selectLSforCutting(std::vector<LevelSetDomain*> &LSDomains)
{
  // select levelsets for cutting
  std::vector<LevelSetDomain*> LSforCutting;
  for (int i=0; i<LSDomains.size(); ++i)
  {
    LevelSetDomain* LS = LSDomains[i];
    if (!LS->meshContainsLS)
      LSforCutting.push_back(LS);
  }
  return LSforCutting;
}

GModel* genLevelSet::cutmesh(GModel* model_, std::vector<LevelSetDomain*> &LSforCutting)
{
  options* opts = new options;
  opts->use_as_library();
  opts->presets(2);

  meshCutter = new MeshCutter(opts);
  meshCutter->load_GModel_mesh(model_);

  // add level-sets
  Mesh* mesh = meshCutter->get_mesh();
  for (int i=0; i<LSforCutting.size(); ++i)
  {
    LevelSetDomain* LS = LSforCutting[i];
    assert (!LS->meshContainsLS);
    addLevelSet(LS);
  }
  // collect all edges before launching process
  mesh->collect_edges();

  meshCutter->process();

  GModel* gm = meshCutter->get_gmodel();
  delete meshCutter;
  meshCutter = NULL;

  return gm;
}

void genLevelSet::addLevelSet(LevelSetDomain* LS)
{
  assert(meshCutter!=NULL);
  Mesh* mesh = meshCutter->get_mesh();
  assert(mesh!=NULL);

  int local_id = mesh->get_entity_manager()->add_face_id(LS->tag);
  mesh->get_entity_manager()->set_real_face(local_id);

  genTerm<genTensor0<double>,0>::ConstHandle D = LS->Distance();

  std::set<MElement*>::const_iterator itep;
  for (itep = LS->pgroup()->begin(); itep != LS->pgroup()->end(); ++itep)
  {
    MElement* ep = *itep;

    int npts = ep->getNumShapeFunctions();
    if (npts==4)
    {
      double u, v, w;
      IntPt pts[npts];
      for (int i=0; i<npts; ++i)
      {
        ep->getNode(i, u, v, w);
        pts[i].pt[0]=u;
        pts[i].pt[1]=v;
        pts[i].pt[2]=w;
        pts[i].weight=1.0;
      }

      std::vector<genScalar<genTensor0<double> > > vvals(npts);
      D->get(ep, npts, &pts[0], vvals);

      int vx[npts];
      double dist[npts];
      for (int i=0; i<npts; ++i)
      {
        MVertex* v = ep->getVertex(i);
        vx[i] = v->getNum();

        dist[i] = vvals[i]();
      }

      LTetra* t = mesh->get_tetra_from_lookup(ep->getNum());
      if (t != NULL)
      {
        for (int i=0; i<npts; ++i)
          assert(t->getVertex(i) == mesh->get_vertex_from_lookup(vx[i]));
        mesh->addLevelset(local_id, t, dist);
      }
    }
  }
}

void genLevelSet::writeCutModel(const std::string &fileName)
{
  if (cutModel!=NULL)
  {
    // write cutModel in MSH file
    char cutModelName[100];
    std::string::size_type idx = fileName.find('.');
    sprintf(cutModelName, "%s_cut.msh", fileName.substr(0, idx).c_str());
    cutModel->writeMSH(cutModelName, 3);
  }
  else
  {
    printf("Impossible to write the cutModel, it does not exist.");
  }
}