// 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>.
//
// Initial design: Frederic Duboeuf (rev.1274)


#include "genEnrichFunctions.h"

//WARNING  GP given in the parametric space of the parent element (cf. MSubElement:getIntegrationPoints)
//         we need to transfert these coordinates
//         in cartesian coordinates using distanceFunction->get

template<class T> void HeavisideFunction<T>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  // We distinguish two cases: internal and external domains
  if (!ele->getParent())
  {
    std::vector<ContainerValType> vvalsd(npts);
    distanceFunction->get(ele,npts,GP,vvalsd);
    for (int i=0; i<npts; ++i)
    {
      if (vvalsd[i]()>0)
        vvals[i]() = 1.;
//    else
//      vvals[i]() = 0.; // default value
    }
  }
  else // Sub-elements are in conformity with the interface
  {
    GaussQuadrature integrator ( GaussQuadrature::Val );
    IntPt* GPEle;
    int nptsEle = integrator.getIntPoints(ele, &GPEle);

    ContainerValType vald;
    distanceFunction->get(ele,nptsEle,GPEle,vald);

    if (vald()>0)
      for (int i=0; i<npts; ++i)
        vvals[i]() = 1.;
//  else
//    for (int i=0; i<npts; ++i)
//      vvals[i]() = 0.; // default value
  }
}


template<class T> void SignFunction<T>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  // We distinguish two cases: internal and external domains
  if (!ele->getParent())
  {
    std::vector<ContainerValType> vvalsd(npts);
    distanceFunction->get(ele,npts,GP,vvalsd);
    for (int i=0; i<npts; ++i)
    {
      if (vvalsd[i]()>0)
        vvals[i]() = 1.;
      else if (vvalsd[i]()<0)
        vvals[i]() =-1.;
//    else
//      vvals[i]() = 0.; // default value
    }
  }
  else // Sub-elements are in conformity with the interface
  {
    GaussQuadrature integrator ( GaussQuadrature::Val );
    IntPt* GPEle;
    int nptsEle = integrator.getIntPoints(ele, &GPEle);

    ContainerValType vald;
    distanceFunction->get(ele,nptsEle,GPEle,vald);

    if (vald()>0)
      for (int i=0; i<npts; ++i)
        vvals[i]() = 1.;
    else if (vald()<0)
      for (int i=0; i<npts; ++i)
        vvals[i]() =-1.;
//  else
//    for (int i=0; i<npts; ++i)
//      vvals[i]() = 0.; // default value
  }
}


template<class T> void GradDiscFunction<T>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  // We need parent parameters
  MElement* elep = ele->getParent();
  if (!elep) elep = ele;

  std::vector<ContainerValType> vvalsd(npts);
  distanceFunction->get(ele,npts,GP,vvalsd);

  double u, v, w;
  int nbSFfct = elep->getNumShapeFunctions();
  IntPt pts[nbSFfct];
  for (int j=0; j<nbSFfct; ++j)
  {
    elep->getNode(j, u, v, w);
    pts[j].pt[0]=u;
    pts[j].pt[1]=v;
    pts[j].pt[2]=w;
    pts[j].weight=1.0;
  }
  std::vector<ContainerValType> vvalsf(nbSFfct);
  distanceFunction->get(elep,nbSFfct,pts,vvalsf);

  for (int i=0; i<npts; ++i)
  {
    u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];
    double val[nbSFfct];
    elep->getShapeFunctions(u, v, w, val);
    for (int j=0; j<nbSFfct; ++j)
    {
      vvals[i]() += ValType(std::fabs(vvalsf[j]()) * val[j]);
    }
    vvals[i]() -= ValType(std::fabs(vvalsd[i]()));                                   // F3 enrichment function suggested in : Moës et al._2003_A computational approach to handle complex microstructure geometries
//  vvals[i]() = (ele->getParent()) ? ValType(std::fabs(vvalsd[i]())) : ValType(1.); // F2 smoothed enrichment function suggested in : Sukumar et al._2001_Modeling holes and inclusions by level sets in the extended finite-element method
//  vvals[i]() = ValType(std::fabs(vvalsd[i]()));                                    // F1 enrichment function suggested in : Belytschko et al._2003_Structured extended finite element methods for solids defined by implicit surfaces
  }
}


template<class T> void GradDiscFunction<T>::getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
{
  // We distinguish two cases: internal and external domains
  std::vector<bool> inExternalDomain(npts,false);
  if (!ele->getParent())
  {
    std::vector<ContainerValType> vvalsd(npts);
    distanceFunction->get(ele,npts,GP,vvalsd);
    for (int i=0; i<npts; ++i)
      if (vvalsd[i]()>0)
        inExternalDomain[i]=true;
  }
  else // Sub-elements are in conformity with the interface
  {
    GaussQuadrature integrator ( GaussQuadrature::Val );
    IntPt* GPEle;
    int nptsEle = integrator.getIntPoints(ele, &GPEle);

    ContainerValType vald;
    distanceFunction->get(ele,nptsEle,GPEle,vald);

    if (vald()>0)
      inExternalDomain = std::vector<bool> (npts,true);
  }

  MElement* elep = ele->getParent();
  if (!elep) elep = ele;

  double u, v, w;
  int nbSFfct = elep->getNumShapeFunctions();
  IntPt pts[nbSFfct];
  for (int j=0; j<nbSFfct; ++j)
  {
    elep->getNode(j, u, v, w);
    pts[j].pt[0]=u;
    pts[j].pt[1]=v;
    pts[j].pt[2]=w;
    pts[j].weight=1.0;
  }
  std::vector<ContainerValType> vvalsf(nbSFfct);
  distanceFunction->get(elep,nbSFfct,pts,vvalsf);

  for (int i=0; i<npts; ++i)
  {
    u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];
    double gradsuvw[nbSFfct][3];
    elep->getGradShapeFunctions(u, v, w, gradsuvw);
    double jac[3][3];
    double invjac[3][3];
    const double detJ = elep->getJacobian(u, v, w, jac);
    inv3x3(jac, invjac);
    if (inExternalDomain[i]) // F3 enrichment function suggested in : Moës et al._2003_A computational approach to handle complex microstructure geometries
    {
      for (int j=0; j<nbSFfct; ++j)
      {
        vgrads[i]() += (ValType(std::fabs(vvalsf[j]())) - vvalsf[j]()) * GradType(
        invjac[0][0] * gradsuvw[j][0] + invjac[0][1] * gradsuvw[j][1] + invjac[0][2] * gradsuvw[j][2],
        invjac[1][0] * gradsuvw[j][0] + invjac[1][1] * gradsuvw[j][1] + invjac[1][2] * gradsuvw[j][2],
        invjac[2][0] * gradsuvw[j][0] + invjac[2][1] * gradsuvw[j][1] + invjac[2][2] * gradsuvw[j][2]);
      }
    }
    else
    {
      for (int j=0; j<nbSFfct; ++j)
      {
        vgrads[i]() += (ValType(std::fabs(vvalsf[j]())) + vvalsf[j]()) * GradType(
        invjac[0][0] * gradsuvw[j][0] + invjac[0][1] * gradsuvw[j][1] + invjac[0][2] * gradsuvw[j][2],
        invjac[1][0] * gradsuvw[j][0] + invjac[1][1] * gradsuvw[j][1] + invjac[1][2] * gradsuvw[j][2],
        invjac[2][0] * gradsuvw[j][0] + invjac[2][1] * gradsuvw[j][1] + invjac[2][2] * gradsuvw[j][2]);
      }
    }
//  if (inExternalDomain[i]) // F2 smoothed enrichment function suggested in : Sukumar et al._2001_Modeling holes and inclusions by level sets in the extended finite-element method
//  {
//    for (int j=0; j<nbSFfct; ++j)
//    {
//      vgrads[i]() += +vvalsf[j]() * GradType(
//      invjac[0][0] * gradsuvw[j][0] + invjac[0][1] * gradsuvw[j][1] + invjac[0][2] * gradsuvw[j][2],
//      invjac[1][0] * gradsuvw[j][0] + invjac[1][1] * gradsuvw[j][1] + invjac[1][2] * gradsuvw[j][2],
//      invjac[2][0] * gradsuvw[j][0] + invjac[2][1] * gradsuvw[j][1] + invjac[2][2] * gradsuvw[j][2]);
//    }
//  }
//  else
//  {
//    for (int j=0; j<nbSFfct; ++j)
//    {
//      vgrads[i]() += -vvalsf[j]() * GradType(
//      invjac[0][0] * gradsuvw[j][0] + invjac[0][1] * gradsuvw[j][1] + invjac[0][2] * gradsuvw[j][2],
//      invjac[1][0] * gradsuvw[j][0] + invjac[1][1] * gradsuvw[j][1] + invjac[1][2] * gradsuvw[j][2],
//      invjac[2][0] * gradsuvw[j][0] + invjac[2][1] * gradsuvw[j][1] + invjac[2][2] * gradsuvw[j][2]);
//    }
//  }
  }
}