// 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>.

#ifndef _MEMBRANE_DOMAIN_H_
#define _MEMBRANE_DOMAIN_H_

#include <string>
#include <vector>

#include "genSupports.h"
#include "genTerm.h"
#include "genTensors.h"
#include "genDataIO.h"
#include "genMatrix.h"
#include "savedGenTerm.h"

// added
#include "membraneStress.h"

class membraneDomain : public genSupport
{
public:
  double _E, _nu; // specific elastic datas (isotropic) shouls not be here (to regroup with the orthotropic case)
  double _thick; // specific to thin bodies
public :
  membraneDomain(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
  }

  template <int n> typename genTerm<genTensor2<double>,n>::Handle Stress(const typename diffTerm<genTensor1<double>,n>::ConstHandle &s);
/*
  template <int n> genTerm<genTensor4<double>,0>::ConstHandle Tangent(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain);
*/
  template <int n1, int n2> typename genTerm<genTensor0<double>, n1+n2>::Handle
  Energy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Strain1, const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain2 );
  membraneDomain() : genSupport() {}
  virtual ~membraneDomain() {};
};

// The Hooke and Tangent tensor can only be build from the element (through a get) GBDTAG
/*
template <int n> genTerm<genTensor4<double>,0>::ConstHandle membraneDomain::Tangent(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain)
{
  return genTerm<genTensor4<double>,0>::ConstHandle(C);
}
*/

template <int n> typename genTerm<genTensor2<double>,n>::Handle membraneDomain::Stress(const typename diffTerm<genTensor1<double>,n>::ConstHandle &sp)
{
  return typename genTerm<genTensor2<double>,n>::Handle(new membraneStress<n>(_E,_nu,_thick,sp));
}


template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
membraneDomain::Energy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Stress_, const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain)
{
  return typename genTerm<genTensor0<double>,n1+n2>::Handle (FullContractedProductCompose(Stress_,Strain)); // no transpose on Strain by construction ??
}

class IsotropicMembraneDomain : public membraneDomain
{
public :
  IsotropicMembraneDomain(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    const double & E_=s.Scalars.find("E")->second;
    const double & nu_=s.Scalars.find("nu")->second;
    const double & h_=s.Scalars.find("thickness")->second;
    _E = E_; _nu = nu_;
    _thick = h_;
  }
  IsotropicMembraneDomain () : membraneDomain() {}
  ~IsotropicMembraneDomain () {}
};
/*
class OrthotropicMembraneDomain : public membraneDomain
{
public :
  std::vector<double> E, nu, G; // E1,E2,E3, nu12,nu13,nu23, G12,G13,G23 (orthotropic_3D)
  genTerm<genTensor4<double>,0>::ConstHandle C; //elastic tensor USED only
public :
  OrthotropicMembraneDomain(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    const std::vector<double> & E_=s.Vectors.find("E")->second;
    const std::vector<double> & nu_=s.Vectors.find("nu")->second;
    const std::vector<double> & G_=s.Vectors.find("G")->second;
    E = E_; nu = nu_; G = G_;
    update();
  }
  OrthotropicMembraneDomain () :  membraneDomain() {}
  ~OrthotropicMembraneDomain () {}
  void update ()
  {
    double S11 = 1. / E[0];
    double S22 = 1. / E[1];
    double S33 = 1. / E[2];
    double S12 = - nu[0] / E[0];
    double S13 = - nu[1] / E[1];
    double S23 = - nu[2] / E[2];

    double dS = S11*S22*S33 - S11*S23*S23 - S22*S13*S13 - S33*S12*S12 + 2*S12*S23*S13;

    genMatrix<double> mC(6,6);
    mC(0,0) = (S22*S33 - S23*S23)/dS;
    mC(1,1) = (S33*S11 - S13*S13)/dS;
    mC(2,2) = (S11*S22 - S12*S12)/dS;
    mC(0,1) = mC(1,0) = (S13*S23 - S12*S33)/dS;
    mC(0,2) = mC(2,0) = (S12*S23 - S13*S22)/dS;
    mC(1,2) = mC(2,1) = (S12*S13 - S23*S11)/dS;
    mC(3,3) = G[2];
    mC(4,4) = G[1];
    mC(5,5) = G[0];

    genTensor4<double> elast;
    FromVoigtNotation(mC,elast);

    C = genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double>>(elast));
  }
};
*/

#endif// _MEMBRANE_DOMAIN_H_