// Bimaterial_xfem - A solver for bimaterial problems using X-FEM
// Copyright (C) 2012-2018 Frederic Duboeuf
//
// See the LICENSE file for license information.
// Please report all bugs and problems to <duboeuf@outlook.com>.

#include <cstring>
#include <cstdio>
#include <fstream>

#include "bimaterialXfemSolver.h"
#include "spaceReducer.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"

#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
#endif


BimaterialXfemSolver::~BimaterialXfemSolver()
{
  for (int i=0; i<LSDomains.size(); ++i) delete LSDomains[i];
  for (int i=0; i<BiDomains.size(); ++i) delete BiDomains[i];
}

void BimaterialXfemSolver::CheckProperties()
{
  ElasticCutmeshXfemSolver::CheckProperties();

  for (int i=0; i<DataLS->LSs.size(); ++i)
  {
    LevelSetDomain* field = new LevelSetDomain(*(DataLS->LSs[i]));
    LSDomains.push_back(field);
  }

  for (unsigned int i = 0; i < DataLS->Domains.size(); ++i)
  {
    if (DataLS->Domains[i]->type=="Bimaterial")
    {
      int dim = DataLS->Domains[i]->dim;
      int tag = DataLS->Domains[i]->tag;
      LevelSetDomain* LS=NULL;

      for (int j=0; j<LSDomains.size(); ++j)
      {
        if ((LSDomains[j]->dim==dim)&&(LSDomains[j]->tag==tag))
          LS=LSDomains[j];
      }

      if (LS)
      {
        BimaterialDomain* field = new BimaterialDomain(*(DataLS->Domains[i]), LS->Distance());
        BiDomains.push_back(field);
      }
      else
        printf("LevelSetDomain with dim=%d and tag=%d is not already existing.\n", dim, tag);
    }
  }
}

void BimaterialXfemSolver::readInputFile ( const std::string &fileName )
{
  modelname = fileName;
  std::cout<<"model name : " << modelname << std::endl ;
  DataLS->readInputFile(fileName);
  CheckProperties();
  printf("--> %d Domains\n", (int)EDomains.size());
  printf("--> %d BCs \n", (int)DataLS->BCs.size());
  printf("--> %d LSs \n", (int)LSDomains.size());
  printf("--> %d Bs \n", (int)BiDomains.size());
}

void BimaterialXfemSolver::CreateFunctionSpaces()
{
  ElasticCutmeshXfemSolver::CreateFunctionSpaces();

  // the dofs of space are numbered first and then enriched.
  // the dofs of enriched space should not be numbered.
  for (int i=0; i<BiDomains.size(); ++i)
  {
    BimaterialDomain* B = BiDomains[i];
    std::cout << "Domain of bi-material interface " << i << " : dim=" << B->dim << ", tag=" << B->tag << ", " << B->group()->psize() << " parent elements" << std::endl;

    std::cout << B->group()->vsize() << " parent vertices" << std::endl;
    std::cout << B->group()->vsize() - B->EnrichedEntities()->size() << " parent vertices along the interface" << std::endl;
    std::cout << B->EnrichedEntities()->size() << " enriched vertices" << std::endl;

    genFilterDofMultiEntities* filterDof = new genFilterDofMultiEntities(B->EnrichedEntities());
    genFSpace<genTensor1<double> >::ConstHandle LagrangeFiltered (new FilteredFSpace<genTensor1<double> >(FSpaceDisp,filterDof));
    genFSpace<genTensor1<double> >::ConstHandle xFemLagrange (new EnrichedFSpace<genTensor1<double> >(LagrangeFiltered, B->EnrichmentFunction()));
    FSpaceDisp = genFSpace<genTensor1<double> >::Handle(new CompositeFSpace<genTensor1<double> >(FSpaceDisp, xFemLagrange));
  }

  genTerm<genTensor2<double>,1>::Handle G(Gradient(FSpaceDisp)*0.5);
  genTerm<genTensor2<double>,1>::Handle CG(new genCache<genTensor2<double>,1>(G));
  EpsilonDisp = (CG+Transpose(CG));
}

PView* BimaterialXfemSolver::buildLevelSetView(const std::string &postFileName, int num)
{
  assert(num<LSDomains.size());
  LevelSetDomain* LS = LSDomains[num];

  genTerm<genTensor0<double>,0>::ConstHandle Field = LS->Distance();

  genGroupOfElements* Support = new genGroupOfElements;
  std::set<MElement*>::iterator itp;
  for (itp=LS->group()->pbegin(); itp!=LS->group()->pend(); ++itp)
  {
    MElement* p = *itp;
    Support->insert(p);
  }
  std::vector<genGroupOfElements*> Supports(1, Support);

  return buildViewNodal(Field, Supports, postFileName);
}

PView* BimaterialXfemSolver::buildEnrichmentFunctionView(const std::string &postFileName, int num)
{
  assert(num<BiDomains.size());
  BimaterialDomain* B = BiDomains[0];

  genTerm<genTensor0<double>,0>::ConstHandle Field = B->EnrichmentFunction();

  std::vector<genGroupOfElements*> Supports;
  Supports.push_back(B->subDomainInt);
  Supports.push_back(B->subDomainExt);

  return buildViewNodal(Field, Supports, postFileName);
}