// 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>.
#ifndef _ELASTICITY_SOLVER_XFEM_H_
#define _ELASTICITY_SOLVER_XFEM_H_

#include <map>
#include <string>
#include "SVector3.h"
#include "dofManager.h"
#include "simpleFunction.h"
#include "functionSpace.h"
#include "gmshLevelset.h"
#include "groupOfElements.h"
#include "mechanicsSolver.h"

#include "mathEvaluator.h"
#include <boost/algorithm/string.hpp>
#include "Composite_Textile/DILevelsetTextile.h"


class function;
class GModel;
class PView;

struct LagrangeMultiplierField {
  int _tag;
  groupOfElements *g;
  double _tau;
  SVector3 _d;
  simpleFunction<double> *_f;
  LagrangeMultiplierField() : g(0),_tag(0){}
};

class LagMultField {

  public :
    int _tag;
    groupOfLagMultElements * _g;
    double _eqfac;
    bool _comp[3];
    simpleFunction<SVector3> * _f;

  private :

  public :

    LagMultField(int tag, int comp, simpleFunction<SVector3> * f, groupOfLagMultElements * g, int eqfac = 1) : _tag(tag), _eqfac(eqfac), _f(f)
    {
      _comp[0] = _comp[1] = _comp[2] = false;
      _comp[comp] =true;
      _g = g;
    }

    void changeLagMultField( simpleFunction<SVector3> * f )
    {
      delete _f;
      _f = f;
    }
};


// Le choix de deux classes dérivées se justifie par la nécéssité de l'appel explicite de SVector3, même avec un template
class MathEvalFunctionSVector3 : public simpleFunction<SVector3>
{
  mutable mathEvaluator _m;
  int _dim;
  public:
  MathEvalFunctionSVector3(std::vector<std::string> &expr,std::vector<std::string> &var) : _m(expr,var), _dim(expr.size()) {}
  virtual SVector3 operator () (double x, double y, double z) const
  { 
    std::vector<double> res(_dim);
    std::vector<double> pos(3);
    pos[0]=x;pos[1]=y;pos[2]=z;
    _m.eval(pos,res);
    if (_dim==3)
      return SVector3(res[0],res[1],res[2]);
    else
      Msg::Error("MathEvalFunction : double operator error (_dim/=3)\n");
  }
  virtual void gradient (double x, double y, double z,
                         SVector3 & dfdx, SVector3 & dfdy, SVector3 & dfdz) const
  { dfdx = dfdy = dfdz = SVector3(.0); }
};

class MathEvalFunctionDouble : public simpleFunction<double>
{
  mutable mathEvaluator _m;
  int _dim;
  public:
  MathEvalFunctionDouble(std::vector<std::string> &expr,std::vector<std::string> &var) : _m(expr,var), _dim(expr.size()) {}
  virtual double operator () (double x, double y, double z) const
  { 
    std::vector<double> res(_dim);
    std::vector<double> pos(3);
    pos[0]=x;pos[1]=y;pos[2]=z;
    _m.eval(pos,res);
    if (_dim==1)
      return double(res[0]);
    else
      Msg::Error("MathEvalFunction : double operator error (_dim/=1)\n");
  }
  virtual void gradient (double x, double y, double z,
                         double & dfdx, double & dfdy, double & dfdz) const
  { dfdx = dfdy = dfdz = double(.0); }
};


// an elastic solver ...
class elasticitySolverXfem : public mechanicsSolver
{
 protected:
      GModel *_pModelUncut;
      // lagrange multipliers function spaces
      FunctionSpace<double> *_LagrangeMultiplierSpace;  // Gaetan method
      std::vector< FunctionSpace < SVector3 > * > _StableLagrangeMultipliersSpace;  // E. Bechet Method
      // simple multiplying function enrichment to enrich the space function related to _lslist
      std::vector<simpleFunctionOnElement<double> *>_funcEnrichmentlist;
      // level set lists to cut the model
      std::vector <gLevelset *> _lslist;
      // level set entity list related to _lslist
      std::vector < std::pair < int,int > > _lsEntitylist;
      // Lagrange multiplier fields
      std::vector<LagrangeMultiplierField> _LagrangeMultiplierFields;  // Gaetan Method
      std::vector<LagMultField> _StableLagrangeMultiplierFields;  // E. Bechet Method
      std::map< int ,std::vector<int> > * _ListHangingNodes;
      // map between uncut and cut physicals
      std::map< int, std::vector< int > > _OriginalPhysicals[4];

 protected :
      // set level sets in the model
      virtual void readLSInInputFile(const std::string &meshFileName);
      // read uncut mesh and cut by level sets _lslist
      virtual void readMeshInInputFile(const std::string &meshFileName);
      // read elastic domain definition and boundary conditions
      virtual void readFieldsInInputFile(const std::string &fn);
      virtual void readBCInInputFile(const std::string &meshFileName);
      // set the map _OriginalPhysicals between uncut and cut physicals
      virtual void setOriginalPhysicals(GModel *pModel, GModel *pModelCut);
      // set the level set vector physicals
      virtual void setLevelSetsPhysicals();
      // build the problem function spaces
      virtual void BuildFunctionSpaces();

 public:
      elasticitySolverXfem(int tag) : mechanicsSolver(tag), _LagrangeMultiplierSpace(0),_ListHangingNodes(0) {}
      virtual void readInputFile(const std::string &meshFileName);
      virtual void setMesh(const std::string &meshFileName);
      virtual void BuildLinearSystem() ;
      virtual PView *buildDisplacementView(const std::string &postFileName);
      virtual PView *buildElasticEnergyView(const std::string &postFileName);
      virtual PView *buildVonMisesView(const std::string &postFileName);
      void addLS(gLevelset * ls){_lslist.push_back(ls);};
      void setUncutGModel(GModel *m){_pModelUncut = m; _dim = 3 ;};
      void addHangingNodes(std::map< int ,std::vector<int> > * ListHangingNodes){_ListHangingNodes = ListHangingNodes;}
      LagMultField* containLagMultField(int physical)
      {
        std::vector<LagMultField>::reverse_iterator itLagMultField = _StableLagrangeMultiplierFields.rbegin();
        for ( ; itLagMultField != _StableLagrangeMultiplierFields.rend(); itLagMultField++)
          if( (*itLagMultField)._tag==physical) return &(*itLagMultField);

        return NULL;
      }
      // std::pair<PView *, PView*> buildErrorEstimateView
      //   (const std::string &errorFileName, double, int);
      // std::pair<PView *, PView*> buildErrorEstimateView
      //   (const std::string &errorFileName, const elasticityData &ref, double, int);
};


#endif