// elastic_genTerm_Cutmesh_Xfem - A linear solver for elastic problems using X-FEM
// Copyright (C) 2010-2026 Eric Bechet, Frederic Duboeuf
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#ifndef _ELASTICXFEM_DOMAIN_H_
#define _ELASTICXFEM_DOMAIN_H_

#include <string>
#include <vector>

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


class StableLagrangeMultiplierDomain : public genSupport
{
public :
  genTensor0<double> eqfac;
  genTensor0<double> k; // stiffness inverse for the condition (0 = infinite stiffness)
  std::set<int> comps;
  genSimpleFunction<genTensor1<double> >* f;

public :
  StableLagrangeMultiplierDomain(const genBoundaryCondition* BC, const genTensor0<double> eqfac_=1.) : eqfac(eqfac_), k(0.), f(0)
  {
    g = new genGroupOfElements(BC->dim, BC->tag);

    dim = BC->dim;
    tag = BC->tag;

    comps.insert(BC->comp1);

    genTensor1<double> fval;
    fval(BC->comp1) = (*BC->fscalar)(0,0,0);
    f = new genSimpleFunctionConstant<genTensor1<double> >(fval);
  }

  void addComponent(const genBoundaryCondition* BC)
  {
    assert(dim == BC->dim);
    assert(tag == BC->tag);

    comps.insert(BC->comp1);

    genTensor1<double> fval = (*f)(0,0,0);
    fval(BC->comp1) = (*BC->fscalar)(0,0,0);
    delete f;
    f = new genSimpleFunctionConstant<genTensor1<double> >(fval);
  }

  template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
  CrossLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &FSpace1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &FSpace2 );

  template <int n1, int n2> typename genTerm<genTensor0<double>, n1+n2>::Handle
  SelfLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &SLag1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &SLag2 );

  template <int n> typename genTerm<genTensor0<double>,n>::Handle Lterm(const typename genTerm<genTensor1<double>,n>::ConstHandle &SLag);

  StableLagrangeMultiplierDomain() : genSupport(), eqfac(1.), k(0.), f(0) {}
  StableLagrangeMultiplierDomain(const StableLagrangeMultiplierDomain &other) : genSupport(other), eqfac(other.eqfac), k(other.k), comps(other.comps), f(new genSimpleFunctionConstant<genTensor1<double> > ((*other.f)(0,0,0))) {}
  virtual ~StableLagrangeMultiplierDomain() {}
};


template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
StableLagrangeMultiplierDomain::CrossLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &FSpace1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &FSpace2)
{
  return typename genTerm<genTensor0<double>,n1+n2>::Handle (Compose<FullContrProdOp>(FSpace1,FSpace2)*eqfac);
}

template <int n1, int n2> typename genTerm<genTensor0<double>,n1+n2>::Handle
StableLagrangeMultiplierDomain::SelfLagTerm( const typename genTerm<genTensor1<double>,n1>::ConstHandle &FSpace1, const typename genTerm<genTensor1<double>,n2>::ConstHandle &FSpace2)
{
  return typename genTerm<genTensor0<double>,n1+n2>::Handle (Compose<FullContrProdOp>(FSpace1,FSpace2)*k); // k peut être un genTerm
}

template <int n> typename genTerm<genTensor0<double>,n>::Handle
StableLagrangeMultiplierDomain::Lterm(const typename genTerm<genTensor1<double>,n>::ConstHandle &FSpaceSLag)
{
  genTerm<genTensor1<double>,0>::Handle Func( new FunctionField<genTensor1<double> >(*f) );
  return typename genTerm<genTensor0<double>,n>::Handle(Compose<FullContrProdOp>(FSpaceSLag,Func)*eqfac);
}


#endif// _ELASTICXFEM_DOMAIN_H_