// Membrane - A linear solver for membrane problem using FEM
// 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>.
// Description: stress for thin bodies formulation

#ifndef _MEMBRANE_STRESS__H_
#define _MEMBRANE_STRESS__H_

#include "genTerm.h"

// Functions call by get with !=n
void stressData(genTensor4<double> &H, genScalar<genTensor2<double> > &Bvals, genScalar<genTensor2<double> > &tmpvals); // n=0
void stressData(genTensor4<double> &H, genVector<genTensor2<double> > &Bvals, genVector<genTensor2<double> > &tmpvals); // n=1


template< int n> class membraneStress : public genTerm< genTensor2<double>, n>
{
public :
  typedef typename ContainerTraits<genTensor2<double>,n>::Container ContainerValType;
  static const int Nsp=n;
protected :
  const double _E;
  const double _nu;
  const double _thick;
  typename diffTerm<genTensor1<double>,n>::ConstHandle _space1;
public :
  membraneStress(const double E, const double nu, const double h, const typename diffTerm<genTensor1<double>,n>::ConstHandle &sp) : _E(E), _nu(nu), _thick(h), _space1(sp){}
  virtual int getNumKeys(MElement* ele,int k=0) const {return _space1->getNumKeys(ele,k);}
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const { _space1->getKeys(ele,keys,k);}
  virtual int getIncidentSpaceTag(int k=0) const { return _space1->getIncidentSpaceTag(k);}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
//  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const {LinearTermBase<typename TensorialTraits<T1>::ValType>::get(ele,npts,GP,vals);}
  virtual membraneStress<n>* clone () const { return new membraneStress<n>(_E,_nu,_thick,_space1);}
  virtual ~membraneStress() {}
};


template<int n> void membraneStress<n>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  vvals.clear(); // Where the values have to be clear ?? GBDTAG

  /* Bvalue */
  std::vector<ContainerValType> Bvals(npts);
  _space1->getgradf(ele,npts,GP,Bvals); // Cache uses and not duplicate the computation ?? GBDTAG


  // stress computation
  for(int j=0;j<npts;j++)
  {
    /* Grad in uvw */ // Hide this in the Functional space GBDTAG ??
    double gradsuvw[256][3];
    ele->getGradShapeFunctions(GP[j].pt[0], GP[j].pt[1], GP[j].pt[2], gradsuvw);

    /* Shell basis vector */ // How to store it and to get access to displacement to compute in the current configuration GBDTAG use the same as in GradMembraneTerm?
    genTensor1<double> phi01(0.,0.,0.);
    genTensor1<double> phi02(0.,0.,0.);
    for(int i=0;i<ele->getNumVertices();i++)
    {
      SVector3 tmp1=(gradsuvw[i][0]*ele->getVertex(i)->point());
      SVector3 tmp2=(gradsuvw[i][1]*ele->getVertex(i)->point());
      genTensor1<double> T1(tmp1[0],tmp1[1],tmp1[2]);
      genTensor1<double> T2(tmp2[0],tmp2[1],tmp2[2]);
      phi01 += T1;
      phi02 += T2;
    }
    // normal
    genTensor1<double> t0 = crossprod(phi01,phi02);
    t0.normalize();
    // dual basis
    genTensor1<double> dp0[2];
    dp0[0] = crossprod(phi02,t0);
    dp0[1] = crossprod(phi01,t0);
    dp0[0] *= (1./dot(dp0[0],phi01));
    dp0[1] *= (1./dot(dp0[1],phi02));

    /* Stiffness */
    genTensor4<double> H(0.); // 3indices but 2 uses GBDTAG ?
    double C11 = _E*_thick/(1.-_nu*_nu);
    double C12 = 0.5*(1.-_nu)*C11;
    C11*=_nu;
    for(int alpha=0;alpha<2;alpha++)
      for(int beta=0;beta<2;beta++)
        for(int gamma=0;gamma<2;gamma++)
          for(int delta=0;delta<2;delta++)
          {
            H(alpha,beta,gamma,delta) += C11*dot(dp0[alpha],dp0[beta])*dot(dp0[gamma],dp0[delta]) + C12*(dot(dp0[alpha],dp0[gamma])*dot(dp0[beta],dp0[delta])+dot(dp0[alpha],dp0[delta])*dot(dp0[beta],dp0[gamma]));
          }
    ContainerValType tmpvals;
    stressData(H,Bvals[j],tmpvals);
    vvals.push_back(tmpvals);
  }
}

#endif //_MEMBRANE_STRESS__H_