// elastic_genTerm - A linear solver for elastic problems 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 _ELASTODYNAMIC_DOMAIN_H_
#define _ELASTODYNAMIC_DOMAIN_H_
#include "genSupports.h"
#include "elasticDomain.h"
#include "genConvert.h"
#include "genCache.h"
#include <cmath>

class ElastodynamicDomain : public genSupport{
public : 
  double omega, rho, E, nu;
  genTerm<genTensor4<std::complex<double> >, 0>::Handle C;
public :
  ElastodynamicDomain( const ElasticDomain &s, const genDomain &p) : genSupport(){
    dim=p.dim;
    tag=p.tag;
    double f = p.Scalars.find("frequence")->second;
    omega = 2*M_PI*f;
    rho = p.Scalars.find("rho")->second;
    C = genTerm<genTensor4<std::complex<double> >, 0>::Handle(new genConvertTerm<genTensor4<std::complex<double> >, genTensor4<double> >(s.C));
  } 
  ElastodynamicDomain( const IsotropicElasticDomain &s, const genDomain &p) : genSupport(){
    dim=p.dim;
    tag=p.tag;
    E = s.E;
    nu = s.nu;
    double f = p.Scalars.find("frequence")->second;
    omega = 2*M_PI*f;
    rho = p.Scalars.find("rho")->second;
    C = genTerm<genTensor4<std::complex<double> >, 0>::Handle( new genConvertTerm<genTensor4<std::complex<double> >, genTensor4<double> >(s.C));
  }
  ElastodynamicDomain( const OrthotropicElasticDomain &s, const genDomain &p) : genSupport(){
    dim=p.dim;
    tag=p.tag;
    double f = p.Scalars.find("frequence")->second;
    omega = 2*M_PI*f;
    rho = p.Scalars.find("rho")->second;
    C = genTerm<genTensor4<std::complex<double> >, 0>::Handle( new genConvertTerm<genTensor4<std::complex<double> >, genTensor4<double> >(s.C));
  }
  template <class GT> typename genTerm<typename GT::ValType, GT::Nsp>::Handle Stress(const std::shared_ptr< GT> &Strain);
  template <int n> genTerm<genTensor4<std::complex<double> >,0>::ConstHandle Tangent(const typename genTerm<genTensor2<std::complex<double> >,n>::ConstHandle &Strain);
  template <class GT1, class GT2> typename genTerm<typename FullContrProdOp< typename GT1::ValType, typename GT2::ValType>::ValType, GT1::Nsp+GT2::Nsp>::Handle
  Energy(const std::shared_ptr< GT1> &Strain1, const std::shared_ptr< GT2> &Strain2);
  ~ElastodynamicDomain() {}
};
  

template <class GT> typename genTerm<typename GT::ValType, GT::Nsp>::Handle ElastodynamicDomain::Stress(const std::shared_ptr<GT> &Strain)
{
  return typename genTerm<typename GT::ValType, GT::Nsp>::Handle(Compose<FullContrProdOp>(C,Strain));
  //return typename genTerm<genTensor2<double>,n>::Handle(FullContractedProductCompose(C,Strain));
}

template <int n> genTerm<genTensor4<std::complex<double> >,0>::ConstHandle ElastodynamicDomain::Tangent(const typename genTerm<genTensor2<std::complex<double> >,n>::ConstHandle &Strain)
{
  return genTerm<genTensor4<std::complex<double> >,0>::ConstHandle(C);
}

template <class GT1, class GT2> typename genTerm<typename FullContrProdOp< typename GT1::ValType, typename GT2::ValType>::ValType, GT1::Nsp + GT2::Nsp>::Handle
ElastodynamicDomain::Energy( const std::shared_ptr<GT1> &Strain1, const std::shared_ptr<GT2> &Strain2)
{
  return typename genTerm<typename FullContrProdOp< typename GT1::ValType, typename GT2::ValType>::ValType ,GT1::Nsp+GT2::Nsp>::Handle (Compose<FullContrProdOp>(Transpose(Strain1),Stress<>(Strain2)));
//  return typename genTerm<double,n1+n2>::Handle (FullContractedProductCompose(Transpose(Strain1),Stress<n2>(Strain2)));
}



#endif // _ELASTODYNAMIC_DOMAIN_H_