// 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: Functional space for membrane problem

#ifndef _MEMBRANE_SPACE_H_
#define _MEMBRANE_SPACE_H_

#include "genFSpace.h"

class membraneSpace : public VectorLagrangeFSpace<double>
{
public:
  typedef TensorialTraits<genTensor1<double> >::ValType ValType;
  typedef TensorialTraits<genTensor1<double> >::GradType GradType;
  typedef TensorialTraits<genTensor1<double> >::HessType HessType;
  typedef ContainerTraits<ValType,1>::Container ContainerValType;
  typedef ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef ContainerTraits<HessType,1>::Container ContainerHessType;

  membraneSpace(const membraneSpace &other) : VectorLagrangeFSpace(other) {}
  membraneSpace() : VectorLagrangeFSpace() {}
  membraneSpace(Along comp1) : VectorLagrangeFSpace(comp1) {}
  membraneSpace(Along comp1, Along comp2) : VectorLagrangeFSpace(comp1,comp2) {}
  membraneSpace(Along comp1, Along comp2, Along comp3) : VectorLagrangeFSpace(comp1,comp2,comp3) {}
  virtual ~membraneSpace() {}

  // overload the gradient operator
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
   // vgrads.clear(); // where the vectors have to be clear GBDTAG??

   // fill gradient value for each Gauss point
    for(int j=0;j<npts;j++)
    {
      /* Grad in uvw */
      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 (ideally should contains shellLocalBasis object) GBDTAG ?
      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;
      }

      /* Gradient value */
      genVector<genTensor2<double> > tmpfull(3*ele->getNumVertices());
      for(int k=0;k<3;k++)
      {
        for(int i=0;i<ele->getNumVertices();i++)
        {
          genTensor2<double> tmptensor(0.);
          tmptensor(0,0) = (phi01[k]*gradsuvw[i][0]);
          tmptensor(0,1) = 0.5*(phi01[k]*gradsuvw[i][1]+phi02[k]*gradsuvw[i][0]);
          tmptensor(1,0) = tmptensor(0,1);
          tmptensor(1,1) = (phi02[k]*gradsuvw[i][1]);
          tmpfull(i+k*ele->getNumVertices()) = tmptensor;
        }
      }
      vgrads[j] =tmpfull;
    }
  }
  virtual membraneSpace* clone () const {return new membraneSpace(*this);}
};
#endif // _MEMBRANE_SPACE_H_