bridging.cpp

Go to the documentation of this file.

/* ./bridging/bridging.cpp 
**********************************
Copyright INRIA and CEA 

author : Guillaume ANCIAUX (anciaux@labri.fr, g.anciaux@laposte.net)

The LibMultiScale is a C++ parallel framework for the multiscale
coupling methods dedicated to material simulations. This framework
provides an API which makes it possible to program coupled simulations
and integration of already existing codes.

This Project is done in a collaboration between INRIA Futurs Bordeaux
within ScAlApplix team and CEA/DPTA Ile de France. 

This software is governed by the CeCILL-C license under French law and
abiding by the rules of distribution of free software.  You can  use, 
modify and/ or redistribute the software under the terms of the CeCILL-C
license as circulated by CEA, CNRS and INRIA at the following URL
"http://www.cecill.info". 

As a counterpart to the access to the source code and  rights to copy,
modify and redistribute granted by the license, users are provided only
with a limited warranty  and the software's author,  the holder of the
economic rights,  and the successive licensors  have only  limited
liability. 

In this respect, the user's attention is drawn to the risks associated
with loading,  using,  modifying and/or developing or reproducing the
software by the user in light of its specific status of free software,
that may mean  that it is complicated to manipulate,  and  that  also
therefore means  that it is reserved for developers  and  experienced
professionals having in-depth computer knowledge. Users are therefore
encouraged to load and test the software's suitability as regards their
requirements in conditions enabling the security of their systems and/or 
data to be ensured and,  more generally, to use and operate it in the 
same conditions as regards security. 

The fact that you are presently reading this means that you have had
knowledge of the CeCILL-C license and that you accept its terms.
***********************************/

//#define TIMER
#define NEED_STAMP
#define NEED_LAMMPS
#define NEED_LIBMESH
#define NEED_MD1D
#define LOCAL_MODULE MOD_COUPLING

#include "../common/common.h"
#include "../dumper/paraview_helper.h"
#include "bridging.h"

#ifdef DEBUG_BRIDGE
unsigned int absolute_lambdasA=0;
unsigned int absolute_lambdasC=0;
unsigned int absolute_positions=0;
#endif

static std::ofstream out_contrainte("contrainte.plot");

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::DumpDistance(){
  if (current_step %100> 0)
    return;

  int taille = atomes_rec.nbElem();
  double xI;
  MechaContLibMesh<Dim> & model = static_cast<MechaContLibMesh<Dim>&>(domC.getModel());  
  _Vec_ & v = domC.getDofs().V();

  std::vector<unsigned int> nodesI;
  std::vector<double> shapesI;

  for (int i = 0 ; i < taille ; ++i){
    xI = atomes_rec.Get(i).position0(0);      
    model.GlobalIndexes(elems_rec.Get(assoc[i]),nodesI);
    model.ComputeShapes(elems_rec.Get(assoc[i]),shapesI,xI);

    out_contrainte << atomes_rec.Get(i).position0(0) << "\t" <<  v(nodesI[0])*shapesI[0] + v(nodesI[1])*shapesI[1] - atomes_rec.Get(i).vitesse(0) << "\t";
  }
  out_contrainte << std::endl;
}


template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::SetParam(char * key,double value){
  if (strcmp(key,"GRID_DIVISIONX")==0){
    division_grille[0] = (int)value;
  }
  if (strcmp(key,"GRID_DIVISIONY")==0){
    division_grille[1] = (int)value;
  }
  if (strcmp(key,"GRID_DIVISIONZ")==0){
    division_grille[2] = (int)value;
  }
  if (strcmp(key,"PARTITION_GENERATION")==0){
    nb_partitions = (int)value;
  }
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::BuildContinuWeight(){

  /* ici demarre la methode de pondération des noeuds du maillage
     qui se trouvent dans la zone de recouvrement.

     On va parcourir les elements si on trouve un element qui
     intersecte la geometrie alors on ajoute celui ci .

     Deuxieme passage pour calculer les pondérations 
     et oui ! on ne connaissait pas la taille avant !

     ZONE : continu 
  */
  
  ContenerElems & c = domC.getContenerElems();
  IteratorElems & it = c.GetIterator();
  RefElt el = it.GetFirst();

  unsigned int nb;

  MechaContLibMesh<Dim> & model = static_cast<MechaContLibMesh<Dim> &>(domC.getModel());
  _Vec_ & p0 = domC.getDofs().P0();

  std::vector< unsigned int > nodes;

  Cube & cube = GeomTools::GetBoundingBox(geom);
  DUMP("la geometrie de la grille ! \n" << cube,DBG_DETAIL);
  grid = new SpatialGrid<int,_Vec_,Dim>(cube,division_grille[0],division_grille[1],division_grille[2]);
  delete &cube;

  Cube & cube1 = GeomTools::GetBoundingBox(geom_fictif);
  DUMP("la geometrie de la grille (fictifs) \n" << cube1 << "\n" << geom_fictif,DBG_DETAIL);
  grid_fictifs = new SpatialGrid<int,_Vec_,Dim>(cube1,division_grille[0],division_grille[1],division_grille[2]);
  delete &cube1;


  //  grid->Dump("grid.vtu");
  //grid_fictifs->Dump("grid-fictifs.vtu");

  int nb_elem = 0;
  int nb_fictifs = 0;

  for(el = it.GetFirst();!it.end();el = it.GetNext())
    {
      model.GlobalIndexes(el,nodes);
      nb = nodes.size();

      bool test,test2;
      test = false;
      test2 = false;

      double x=0.0,y=0.0,z=0.0;
     
      for (unsigned int i=0; i< nb ; ++i)
        {
          DUMP("invalid mass ? " << nodes[i] << " " <<  model.Mass(nodes[i]),DBG_ALL);

          if(model.Mass(nodes[i]) == 0)
            continue;

          x = p0[nodes[i]*Dim];
          if (Dim > 1)
            y = p0[nodes[i]*Dim+1];
          if(Dim == 3)
            z = p0[nodes[i]*Dim+2];
          
          test = test | geom.Contient(x,y,z);
          test2 = test2 | geom_fictif.Contient(x,y,z);

          DUMP("is " << x << " " << y << " " << z << " inside " << geom_fictif << " ? result = " << test2,DBG_ALL);
        }

      if (Dim == 3)
        if (nb < 4){
          test = false;
          test2 = false;
        }

      if(test){
        DUMP("element(" << nb_elem << ") added nodes = [ " << nodes[0] << " " << nodes[1] << " " << nodes[2],DBG_ALL);
        STARTTIMER("Filling spatial-grid");     
        grid->add_finiteElement(nb_elem,p0,nodes);
        STOPTIMER("Filling spatial-grid");      
        //ajout concret du noeud dans la liste des noeuds
        for (unsigned int i=0;i < nb;++i){
          x = p0[nodes[i]*Dim];
          if (Dim > 1)
            y = p0[nodes[i]*Dim+1];
          if(Dim == 3)
            z = p0[nodes[i]*Dim+2];
          
          //j'agrandi la geometrie locale
          LocalGeom.ExtendBoundingBox(x,y,z);

          if (noeuds_rec.Recherche(nodes[i])==-1)
            noeuds_rec.Ajouter(nodes[i]);
          
          DUMP("Ajoute le noeuds a la position " << x << " " << y << " " << z,DBG_ALL);
        }
        ++nb_elem;
        elems_rec.Ajouter(el);
      }
      
      if (test2){
        grid_fictifs->add_finiteElement(nb_fictifs,p0,nodes);
        for (unsigned int i=0;i < nb;++i){
          x = p0[nodes[i]*Dim];
          if (Dim > 1)
            y = p0[nodes[i]*Dim+1];
          if(Dim == 3)
            z = p0[nodes[i]*Dim+2];


          if (noeuds_fictifs.Recherche(nodes[i])==-1)
            noeuds_fictifs.Ajouter(nodes[i]);


          DUMP("Ajoute le noeuds fictif a la position " << x << " " << y << " " << z,DBG_ALL);
        }
        ++nb_fictifs;
        elems_fictifs.Ajouter(el);
      }
    }
      
  delete &it;
  
  
  DUMP("on a trouve " << nb_fictifs << " elements concernes(fictifs)",DBG_INFO);
  DUMP("on a en moyenne " << grid->getAverageEltByBlock() << " elements par block",DBG_INFO);
  DUMP("on a en moyenne " << grid_fictifs->getAverageEltByBlock() << " elements par block(fictifs)",DBG_INFO);
  // grid->Dump("grid.vtu");
  //grid_fictifs->Dump("grid-fictifs.vtu");
  nb = noeuds_rec.nbElem();
    DUMP("on a trouve " << nb_elem << " elements concernes (nb = " << nb << ")",DBG_INFO);
  if (!nb){DUMP("On a trouve aucun noeud dans le recouvrement",DBG_INFO);return;}
  DUMP("On a trouve " << nb << " noeuds concernes dans " << elems_rec.nbElem() << " elements",DBG_INFO);
  /* j'alloue ici le vecteur alpha pour continu */


  lambdasC = new double[nb];
#ifdef DEBUG_BRIDGE
  absolute_lambdasC = nb;
#endif


  IteratorNodesRec & itrec = noeuds_rec.GetIterator();
  RefNode n = itrec.GetFirst();
  unsigned int j = 0;

  while (!itrec.end())
    {
#ifdef DEBUG_BRIDGE
      if (j >= absolute_lambdasC) FATAL("depassement detecte");
#endif
      if (Dim == 1)
        lambdasC[j] = P.weight(CONTINUFLAG,p0[n]);
      if (Dim == 2)
        lambdasC[j] = P.weight(CONTINUFLAG,p0[n*Dim],p0[n*Dim+1]);
      if (Dim == 3)
        lambdasC[j] = P.weight(CONTINUFLAG,p0[n*Dim],p0[n*Dim+1],p0[n*Dim+2]);

      DUMP("node " << n << " placed in position " << j,DBG_ALL);

      n = itrec.GetNext();
      ++j;
    }


  
  delete &itrec;
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::BuildAtomsWeight(){
  /* ici demmarre la fonction de pondération des atomes.
     on parcours de la meme manière les atomes .
     s'ils sont dans la geometrie on les notes

     deuxieme passage pour ponderer

     ZONE : atoms

  */

  DUMP("la geometrie du rec (fictifs) \n" << geom_fictif,DBG_INFO);

  DUMP("taille de atomes_rec avant de demarrer " << atomes_rec.nbElem(),DBG_INFO);

  ContenerAtoms & c = domA.getContener();
  IteratorAtoms& it = c.GetIterator();
  RefAtom at = it.GetFirst();
  //  int index = 0;
  bool test,test2;

  while (!it.end())
    {
      DUMP("Testing atom at position " << FORMATREAL(at.position0(0)),DBG_ALL);
      
      test = false;
      test2 = false;
      
      if (Dim == 1){
        if (geom.Contient(at.position0(0)))
          test = true;
        if (geom_fictif.Contient(at.position0(0)))
          test2 = true;
      }

      if (Dim == 2){
        if (geom.Contient(at.position0(0),at.position0(1)))
          test = true;
        if (geom_fictif.Contient(at.position0(0),at.position0(1)))
          test2 = true;
      }

      
      if (Dim == 3){
        if (geom.Contient(at.position0(0),at.position0(1),at.position0(2)))
          test = true;
        if (geom_fictif.Contient(at.position0(0),at.position0(1),at.position0(2)))
          test2 = true;
      } 

      if (test){
        DUMP("Accepting atom at position " << FORMATREAL(at.position0(0)),DBG_ALL);
        atomes_rec.Ajouter(at);
        //j'agrandi la geometrie locale
        if(Dim == 1)
          LocalGeom.ExtendBoundingBox(at.position0(0),0,0);
        if(Dim == 2)
          LocalGeom.ExtendBoundingBox(at.position0(0),at.position0(1),0);
        if(Dim == 3)
          LocalGeom.ExtendBoundingBox(at.position0(0),at.position0(1),at.position0(2));
        //      ++index;
      }
      if (test2){
        atomes_fictifs.Ajouter(at);
        if(Dim == 1)
          LocalGeom.ExtendBoundingBox(at.position0(0),0,0);
        if(Dim == 2)
          LocalGeom.ExtendBoundingBox(at.position0(0),at.position0(1),0);
        if(Dim == 3)
          LocalGeom.ExtendBoundingBox(at.position0(0),at.position0(1),at.position0(2));
        DUMP("Accepting atom at position " << at.position0(0) << " " << at.position0(1) << " in boundary",DBG_ALL);  
      }
      
      at = it.GetNext();
    }

  int nb = atomes_rec.nbElem();
  if (!nb){DUMP("On a trouve aucun atome dans le recouvrement",DBG_INFO);return;}
  DUMP("On a touve " << nb << " atomes concernes et " <<  atomes_fictifs.nbElem() << " atomes qu'il faudra corriger de l'effet de surface",DBG_INFO);

  /* j'alloue ici le vecteur de poids */
  lambdasA = new double[nb+PADDING_ATOMS];
#ifdef DEBUG_BRIDGE
  absolute_lambdasA = nb+PADDING_ATOMS;
#endif
  IteratorAtomsRec & itrec = atomes_rec.GetIterator();
  at = itrec.GetFirst();
  unsigned int i = 0;

  while (!itrec.end())
    {
#ifdef DEBUG_BRIDGE
      if (i >= absolute_lambdasA)
        FATAL("depassement detecte");
#endif

      DUMP("ponderation of rec atom " << i,DBG_ALL);
      if (Dim == 1)
        lambdasA[i] = P.weight(ATOMEFLAG,at.position0(0));
      if (Dim == 2)
        lambdasA[i] = P.weight(ATOMEFLAG,at.position0(0),at.position0(1));
      if (Dim == 3)
        lambdasA[i] = P.weight(ATOMEFLAG,at.position0(0),at.position0(1),at.position0(2));

      at = itrec.GetNext();
      ++i;
    }

  delete &itrec;
  delete &it;
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::BuildPositions(){

  if (atomes_rec.nbElem() == 0){
    positions = NULL;
    return;
  }
  positions = new double[Dim*atomes_rec.nbElem()];
  DUMP("making positions for " << atomes_rec.nbElem() << " atoms in rec",DBG_INFO);
#ifdef DEBUG_BRIDGE
  absolute_positions = atomes_rec.nbElem();
  unsigned int nmax = atomes_rec.nbElem();
#endif
  for (unsigned int i = 0 ; i < atomes_rec.nbElem() ; ++i)
    {
#ifdef DEBUG_BRIDGE
      if (i >= nmax) FATAL("depassement detecte");
#endif
      positions[Dim*i] =  atomes_rec.Get(i).position0(0);
      if (Dim > 1)
        positions[Dim*i + 1] = atomes_rec.Get(i).position0(1);
      if (Dim == 3)
        positions[Dim*i + 2]= atomes_rec.Get(i).position0(2);
    }
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::BuildFictifsPositions(){
  if (atomes_fictifs.nbElem() == 0){
    positions = NULL;
    return;
  }
  positions = new double[Dim*atomes_fictifs.nbElem()];
  DUMP("making positions for " << atomes_fictifs.nbElem() << " atoms in fictifs",DBG_INFO);
#ifdef DEBUG_BRIDGE
  unsigned int nmax = atomes_fictifs.nbElem();
  absolute_positions = atomes_fictifs.nbElem();
#endif

  for (unsigned int i = 0 ; i < atomes_fictifs.nbElem() ; ++i)
    {
#ifdef DEBUG_BRIDGE
      if (i >= nmax) FATAL("depassement detecte");
#endif
      positions[Dim*i] =  atomes_fictifs.Get(i).position0(0);
      if (Dim > 1)
        positions[Dim*i + 1] = atomes_fictifs.Get(i).position0(1);
      if (Dim == 3)
        positions[Dim*i + 2]= atomes_fictifs.Get(i).position0(2);
    }
}


template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::BuildShapeMatrix(double * positions,unsigned int nb_atoms){

  /* ici demarre l'algo de recoupement des éléments 

  On associe dans un tableau de taille nbElementRec 
  des petites matrices de taille nb_atome_par_element x connectivite 
  
  qui vont conenir un bloc de la  ShapeMatrix 

  pour cela on commence par creer un tableau qui donne l'element conteneur 
  pour chaque atome

  ensuite il nous faut le nombre de d'atomes par elements

  ensuite on parcours les elements et on creer la matrice
  aux bonnes dimensions et on la remplie !!

  ZONE : FE
  */

  std::vector<double> atomes_probleme;

  MechaContLibMesh<Dim> & model = static_cast<MechaContLibMesh<Dim>&>(domC.getModel());

  assoc.resize(nb_atoms);
  memset(&assoc[0],-1,sizeof(int)*nb_atoms);

  if (elems_rec.nbElem() == 0){
    DUMP("elem_rec est vide",DBG_WARNING);
    return;
  } 

  //  unsigned int j=0;
  
  std::vector<unsigned int> nb_atome_par_element;
  nb_atome_par_element.resize(elems_rec.nbElem());  

  /* construction du tableau d'assoc atoms elem */

  bool test;

  STARTTIMER("Construction association");

  for (unsigned int i = 0 ; i < nb_atoms ; ++i){

#ifndef TIMER
    if(i%100000==0){
      DUMP("construction de l'association - atome " << i << "/" << nb_atoms,DBG_INFO);
    }
#endif

#ifdef DEBUG_BRIDGE
    if (i >= absolute_positions) FATAL("depassement detecte");
#endif
    
    double x = positions[i*Dim];
    double y = 0.0;
    double z = 0.0;
    if (Dim > 1)
      y = positions[i*Dim+1];
    if (Dim == 3)
      z = positions[i*Dim+2];

    std::set<int> & subset_elts = grid->findSet(x,y,z);
    std::set<int>::iterator it = subset_elts.begin(); 

    for (it = subset_elts.begin(); it != subset_elts.end();++it){
      
      typename Bridging<DomainA,DomainC,Dim,Pond>::RefElt el = elems_rec.Get(*it);
      unsigned int j=*it;

      test = false;
      
      if (Dim == 1)
        if (el.Contient(x))
          test = true;
      
      if (Dim == 2)
        if (el.Contient(x,y))
          test = true;
      
      if (Dim == 3)
        if (el.Contient(x,y,z))
          test = true;
      
      if (test)
        {
          if (assoc[i]!= -1){FATAL("gros soucy qd meme ca fait deux fois que je trouve le meme atome");}
        
          assoc[i] = j;

//        if (i == 39) {
//          DUMP("moi j'ai l'atome 39 (x = " << x << " y = " << y << " z = " << z << ") ! et voici l'element ");
//          MechaContLibMesh<Dim> & model = static_cast<MechaContLibMesh<Dim> &>(domC.getModel());
//          _Vec_ & p0 = domC.getDofs().P0();
            
//          std::vector< unsigned int > nodes;

//          model.GlobalIndexes(el,nodes);
//          unsigned int nb = nodes.size();

//          double x=0.0,y=0.0,z=0.0;
     
//          for (unsigned int i=0; i< nb ; ++i)
//            {
//              x = p0[nodes[i]*Dim];
//              y = p0[nodes[i]*Dim+1];
//              z = p0[nodes[i]*Dim+2];
//              DUMP("node "  << i << " x = " << x << " y = " << y 
//                   << " z = " << z);
//            }
//        }
          ++assoc_found;
          break;
        }

    }
    
    if (it == subset_elts.end()){
//       int a=1;
//       while (a){}
      //      DumpAtome(x,y,z);
      //FATAL("l atome " << i << " du recouvrement a la position " << x 
        //    << "," << y << "," << z << " n'a pas trouve son element");

      assoc[i]=-1;
    }
  }
  STOPTIMER("Construction association");

  delete grid;

  /* construction du nombre d'atomes par elements */ 

  for (unsigned int i = 0 ; i < nb_atoms ; ++i){
    if (assoc[i]!=-1)
      ++nb_atome_par_element[assoc[i]];
  }  

  /* allocation des subshapematrices */

  filterContenerElemsAndNodes(nb_atome_par_element,elems_rec,noeuds_rec,assoc,lambdasC);

  //  for (unsigned int i=0;i < elems_rec.nbElem() ; ++i){std::cout << "nb_atome_par_element["<< i << "] = " << nb_atome_par_element[i] << endl;}

  typename Bridging<DomainA,DomainC,Dim,Pond>::IteratorElemsRec & itc = elems_rec.GetIterator();
  typename Bridging<DomainA,DomainC,Dim,Pond>::RefElt el = itc.GetFirst();

  Smatrix = new ShapeMatrix<_Vec_,Dim>(elems_rec.nbElem());
  el = itc.GetFirst();
  unsigned int i = 0;
  while (!itc.end()){
    DUMP("declare la smatrice numero " << i << " qui sera de taille " << nb_atome_par_element[i] << "x" << el.n_nodes(),DBG_ALL);
    if (nb_atome_par_element[i]>0)
      Smatrix->SetSub(i,nb_atome_par_element[i],el.n_nodes());
    ++i;
    el = itc.GetNext();
  }  

  //fait l'allocation effective qui sera contigue (de preference !!)
  Smatrix->Allocate();

  /* mise en place de l'indirection et les shapes*/
  std::vector<double> shapes;
  std::vector<unsigned int> nodes;
  int indirection=0;

  for (unsigned int i = 0 ; i < nb_atoms ; ++i){
    if (assoc[i]==-1)
      continue;
    
    
//     if (i%10000 == 0)
//       DUMP("computing shapes " << i << "/" << nb_atoms);

    double x = positions[i*Dim];
    double y = 0.0;
    double z = 0.0;

    if (Dim > 1)
      y = positions[i*Dim+1];
    if (Dim == 3)
      z = positions[i*Dim+2];


    el = elems_rec.Get(assoc[i]);

    if (Dim == 1)
      model.ComputeShapes(el,shapes,x);

    if (Dim == 2)
      model.ComputeShapes(el,shapes,x,y);

    if (Dim == 3)
      model.ComputeShapes(el,shapes,x,y,z);

    model.GlobalIndexes(el,nodes);

    for(unsigned int j=0;j<shapes.size();++j){
      DUMP("inserting shape of node " << nodes[j] << "(part of elem " << assoc[i] << ") for atom " << i << " (index " << indirection << ")",DBG_ALL);
      (*Smatrix)(indirection,nodes[j],noeuds_rec.Recherche(nodes[j]),assoc[i],shapes[j]);
    }
    ++indirection;
  }  
  
  delete &itc;
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::BuildShapeMatrixFictifs(double * positions,unsigned int nb_atoms){
  /* meme chose mais pour les atomes fictifs seulement 
  ZONE : FE
  */

  MechaContLibMesh<Dim> & model = static_cast<MechaContLibMesh<Dim>&>(domC.getModel());

 
  assoc_fictifs.resize(nb_atoms);
  memset(&assoc_fictifs[0],-1,sizeof(int)*nb_atoms);

  if (elems_fictifs.nbElem() == 0){
    DUMP("elem_fictifs est vide",DBG_WARNING);
    return;
  } 

  std::vector<unsigned int> nb_atome_par_element;
  nb_atome_par_element.resize(elems_fictifs.nbElem());  

  /* construction du tableau d'assoc atoms elem */

  bool test;

  for (unsigned int i = 0 ; i < nb_atoms ; ++i){

    if(i%10000==0){
      DUMP("construction de l'association - atome(fictifs) " << i << "/" << nb_atoms,DBG_INFO);
    }


    double x = positions[i*Dim];
    double y = 0.0;
    double z = 0.0;
    if (Dim > 1)
      y = positions[i*Dim+1];
    if (Dim == 3)
      z = positions[i*Dim+2];

    std::set<int>::iterator it;
    std::set<int> & subset_elts = grid_fictifs->findSet(x,y,z);
    for (it = subset_elts.begin(); it != subset_elts.end();it++){

      DUMP("accessing element " << *it << " in subset",DBG_ALL);
      typename Bridging<DomainA,DomainC,Dim,Pond>::RefElt el = elems_fictifs.Get(*it);
      unsigned int j=*it;

      test = false;

      if (Dim == 1)
        if (el.Contient(x))
          test = true;

      if (Dim == 2)
        if (el.Contient(x,y))
          test = true;

      if (Dim == 3)
        if (el.Contient(x,y,z))
          test = true;

      if (test)
        {
          assoc_fictifs[i] = j;
          ++assoc_found_fictifs;
          break;
        }
    }
    if (it == subset_elts.end()){
      //      DumpAtome(x,y,z);
      //FATAL("l'atome " << i << " du recouvrement a la position " << x << " " << y << " " << z
        //    << " n'a pas trouve son element");
      assoc_fictifs[i] = -1;
    }
    
    
  }  

  delete grid_fictifs;
  
  /* construction du nombre d'atomes par elements */ 
  
  for (unsigned int i = 0 ; i < nb_atoms ; ++i){
    if (assoc_fictifs[i]!=-1)
      ++nb_atome_par_element[assoc_fictifs[i]];
  }  

  filterContenerElemsAndNodes(nb_atome_par_element,elems_fictifs,noeuds_fictifs,assoc_fictifs,NULL);  

  /* allocation des subshapematrices */

  typename Bridging<DomainA,DomainC,Dim,Pond>::IteratorElemsRec & itc = elems_fictifs.GetIterator();
  typename Bridging<DomainA,DomainC,Dim,Pond>::RefElt el = itc.GetFirst();

  Smatrix_fictifs = new ShapeMatrix<_Vec_,Dim>(elems_fictifs.nbElem());
  el = itc.GetFirst();
  int i = 0;
  while (!itc.end()){
    Smatrix_fictifs->SetSub(i,nb_atome_par_element[i],el.n_nodes());
    DUMP("declare la smatrice numero " << i << " qui sera de taille " << nb_atome_par_element[i] << "x" << el.n_nodes(),DBG_DETAIL);
    ++i;
    el = itc.GetNext();
  }  

  //fait lallocation reelle des coeficients
  Smatrix_fictifs->Allocate();

  /* mise en place de l'indirection et les shapes*/
  std::vector<double> shapes;
  std::vector<unsigned int> nodes;
  int indirection=0;

  for (unsigned int i = 0 ; i < nb_atoms ; ++i){
    if(assoc_fictifs[i]==-1)
      continue;

//     if (i%10000 == 0)
//       DUMP("computing shapes " << i << "/" << nb_atoms);

    double x = positions[i*Dim];
    double y = 0.0;
    double z = 0.0;

    if (Dim > 1)
      y = positions[i*Dim+1];
    if (Dim == 3)
      z = positions[i*Dim+2];

    
    el = elems_fictifs.Get(assoc_fictifs[i]);

    if (Dim == 1)
      model.ComputeShapes(el,shapes,x);

    if (Dim == 2)
      model.ComputeShapes(el,shapes,x,y);

    if (Dim == 3)
      model.ComputeShapes(el,shapes,x,y,z);

    model.GlobalIndexes(el,nodes);

    for(unsigned int j=0;j<shapes.size();++j){
      DUMP("inserting shape of node " << nodes[j] << "(part of elem " << assoc_fictifs[i] << ") for atom " << i,DBG_ALL);
      (*Smatrix_fictifs)(indirection,nodes[j],noeuds_fictifs.Recherche(nodes[j]),assoc_fictifs[i],shapes[j]);
    }
    ++indirection;
  }
  
  delete &itc;
}


template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::CorrectSurfaceEffect(){

  if (elems_fictifs.nbElem() == 0){
    DUMP("Warning : elems rec contener is empty -> atoms surface effect correction will not be made(check geometries)",DBG_WARNING);
    return;
  }

  if (atomes_fictifs.nbElem() == 0){
    DUMP("Warning : atoms rec contener is empty -> atoms surface effect correction will not be made(check geometries)",DBG_WARNING);
    return;
  }


  /* Dans la geometrie Boule je met les atomes soumis aux effets de surface */
  /* je les interpole du maillage justement pour eviter les effets indésirables */

  _Vec_ & v = domC.getDofs().V();
  _Vec_ & u = domC.getDofs().U();
  _Vec_ & p0 = domC.getDofs().P0();
  
  typename Bridging<DomainA,DomainC,Dim,Pond>::IteratorAtomsRec & itfictif = atomes_fictifs.GetIterator();
  typename Bridging<DomainA,DomainC,Dim,Pond>::RefAtom at = itfictif.GetFirst();
  unsigned int at_index = 0;

  int * assoc_tmp = &assoc_fictifs[0];

  for (at = itfictif.GetFirst(); !itfictif.end() ; ++at_index, at = itfictif.GetNext())
    {
      //je zappe un atome a indirection invalide
      if (assoc_tmp[at_index] == -1){
        DUMP("Warning : dofs conteners of constraint system have not been filtered (index = " << at_index << ")",DBG_WARNING);
        continue;
      }
        
      /* je fait du dump */
      at.force(0) = 0;
      at.vitesse(0) = Smatrix_fictifs->InterpolateAtom(at_index,v,assoc_tmp[at_index],0);
      at.position(0) = Smatrix_fictifs->InterpolateAtom(at_index,p0,assoc_tmp[at_index],0) + 
        // at.position0(0) +
        Smatrix_fictifs->InterpolateAtom(at_index,u,assoc_tmp[at_index],0);

      DUMP("correcting surface atom " << at_index 
           << " for Y dimension with " << Smatrix_fictifs->InterpolateAtom(at_index,v,assoc_tmp[at_index],1),DBG_ALL);
      
      if (Dim > 1){
        at.force(1) = 0;
        at.vitesse(1) = Smatrix_fictifs->InterpolateAtom(at_index,v,assoc_tmp[at_index],1);
//      at.position(1) = //Smatrix_fictifs->InterpolateAtom(at_index,p0,assoc_tmp[at_index],1)
//        at.position0(1) + Smatrix_fictifs->InterpolateAtom(at_index,u,assoc_tmp[at_index],1);
      }//  else{
//      at.force(1) = 0;
//      at.vitesse(1) = 0;
//      at.position(1) = 0;
//       }
      
      if (Dim == 3){
        at.force(2) = 0;
        at.vitesse(2) = Smatrix_fictifs->InterpolateAtom(at_index,v,assoc_tmp[at_index],2);
//      at.position(2) = //Smatrix_fictifs->InterpolateAtom(at_index,p0,assoc_tmp[at_index],2) 
//        at.position0(2) + Smatrix_fictifs->InterpolateAtom(at_index,u,assoc_tmp[at_index],2);
      }
//       else{
//      at.force(2) = 0;
//      at.vitesse(2) = 0;
//      at.position(2) = 0;
//       }
      
//       at = itfictif.GetNext();
//       ++at_index;
    }

  delete &itfictif;
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::filterContenerElemsAndNodes(std::vector<unsigned int> & nb_atome_par_element,
                                                                   ContenerElemsRec & elems,
                                                                   ContenerNodesRec & nodes,
                                                                   std::vector<int> & assoc,double * lambda){

  //  MechaContLibMesh<Dim> & model = static_cast<MechaContLibMesh<Dim> &>(domC.getModel());

  ContenerElemsRec new_t;
  ContenerNodesRec new_t_nodes;
  std::vector<int> new_index;
  std::vector<unsigned int> nodeIndexes;

  lambda = NULL;
        
  for (unsigned int i = 0 ; i <  nb_atome_par_element.size(); ++i){
    DUMP("filtering odl index element " << i,DBG_ALL);
    if (nb_atome_par_element[i] > 0){
      DUMP("for contener, elem " << i << " becomes " << new_t.nbElem(),DBG_ALL);
      new_index.push_back(new_t.nbElem());

      DUMP("maintenant nb_atome_par_element[" << new_t.nbElem() << "]= " << nb_atome_par_element[i] << " (i=" << i << ")",DBG_ALL);
      nb_atome_par_element[new_t.nbElem()] = nb_atome_par_element[i];
      new_t.Ajouter(elems.Get(i));

//       RefElt elem = elems.Get(i);
//       model.GlobalIndexes(elem,nodeIndexes);

//       for (unsigned int j = 0 ; j < nodeIndexes.size() ; ++j){
//      vector<unsigned int>::iterator it;
//      if(new_t_nodes.Recherche(nodeIndexes[j])==-1){
//        if (lambda)
//          lambda[new_t_nodes.nbElem()] = lambda[nodes.Recherche(nodeIndexes[j])];

//        new_t_nodes.Ajouter(nodeIndexes[i]);
//      }
//       }

    }
    else
      new_index.push_back(new_t.nbElem());
  }

  nb_atome_par_element.resize(new_t.nbElem());
  
  elems.Empty();
  //  nodes.Empty();
  
  for (unsigned int i = 0 ; i < new_t.nbElem();++i)
    elems.Ajouter(new_t.Get(i));

//   for (unsigned int i = 0 ; i < new_t_nodes.nbElem();++i)
//     nodes.Ajouter(new_t_nodes.Get(i));
    
  for (unsigned int i = 0 ; i < assoc.size();++i){
    if (assoc[i] != -1)
      assoc[i] = new_index[assoc[i]];
  }
  

}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::filterContener(ContenerAtomsRec & atoms,std::vector<int> & v){
  //TODO fair un algo un peu mieux que cette merde de copie
  ContenerTableau<RefAtom> new_t;

  for (unsigned int i = 0 ; i < v.size();++i)
    if(v[i]!=-1){
      DUMP("for contener, atom " << i << " becomes " << new_t.nbElem(),DBG_ALL);
      new_t.Ajouter(atoms.Get(i));
    }

  atoms.Empty();
  for (unsigned int i = 0 ; i < new_t.nbElem();++i)
    atoms.Ajouter(new_t.Get(i));

}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
#ifdef DEBUG_BRIDGE
void Bridging<DomainA,DomainC,Dim,Pond>::filterArray(double * array,unsigned int nmax,std::vector<int> & v){
#else
void Bridging<DomainA,DomainC,Dim,Pond>::filterArray(double * array,std::vector<int> & v){
#endif

  unsigned int index=0;

  for (unsigned int i = 0 ; i < v.size();++i)
    if(v[i]!=-1){
      DUMP("for contener, atom " << i << " becomes " << index,DBG_ALL);
#ifdef DEBUG_BRIDGE
      if (i >= nmax || index >= nmax)
        FATAL("depassement detecte");
#endif
      array[index] = array[i];
      ++index;
    }
}


template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::rebuildAssoc(std::vector<int> & v){
  //TODO fair un algo un peu mieux que cette merde de copie
  std::vector<int> tmp;
  for (unsigned int i = 0 ; i < v.size();++i)
    if(v[i]!=-1){
      DUMP("for assoc, atom " << i << " becomes " << tmp.size() << " value = (" << v[i] << ")",DBG_ALL);
      tmp.push_back(v[i]);
    }

  v.clear();
  for (unsigned int i = 0 ; i < tmp.size();++i){
    DUMP("pushing " << tmp[i] << " at position " << i,DBG_ALL);    
    v.push_back(tmp[i]);
  }

}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::Dump(char * prefix){
  char filename[256];

  snprintf(filename,256,"%s-%s-%s",name,prefix,"recatomes");
  Dump(filename,atomes_rec,true);
  snprintf(filename,256,"%s-%s-%s",name,prefix,"recelems");
  Dump(filename,noeuds_rec,elems_rec,true);
  snprintf(filename,256,"%s-%s-%s",name,prefix,"fictifsatomes");
  Dump(filename,atomes_fictifs,false);
  snprintf(filename,256,"%s-%s-%s",name,prefix,"fictifselems");
  Dump(filename,noeuds_fictifs,elems_fictifs,false);
}



template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::Dump(char * fichier,ContenerTableau<RefAtom> & atomes,bool dumpWeights){

  UnitsConverter unit;
  unit.SetReadUnits(UnitsConverter::AtomsUnits);


  /***************************************************/
  /* je fait le dump des poids des atomes */
  /***************************************************/

  //  IteratorAtomsRec & itrec = atomes_fictifs.GetIterator();
  if (atomes.nbElem() == 0){
    DUMP("Warning : not dumping to file Bridged information as atoms contener seems to be empty",DBG_WARNING);
    return;
  }

  IteratorAtomsRec & itrec = atomes.GetIterator();

  itrec.GetFirst();
  RefAtom at = itrec.GetFirst();

  LMFile file;
  unsigned int nb = atomes.nbElem();
  //unsigned int nb = atomes_fictifs.nbElem();
  
  char filename[256];
  sprintf(filename,"%s-atoms_weigth.vtu",fichier);
  file.open(filename);

  ParaviewHelper paraHelper;
  paraHelper.setOutputFile(file);
  paraHelper.SetMode(TEXT);

  paraHelper.write_header(nb,1);
  paraHelper.startDofList();
 

  while (!itrec.end())
    {
      double x = 0, y = 0 , z = 0;
      
      x = at.position(0)*unit.etalon_metre;
      if (Dim > 1)
        y = at.position(1)*unit.etalon_metre;
      if (Dim > 2)
        z = at.position(2)*unit.etalon_metre;

      paraHelper.pushDouble(x);
      paraHelper.pushDouble(y);
      paraHelper.pushDouble(z);
      at = itrec.GetNext();
    }

  paraHelper.endDofList();
  paraHelper.startCellsConnectivityList();
  
  for (unsigned int i=0;i < nb;++i)
    {
      paraHelper.pushInteger(i);      
    }    
  
  paraHelper.endCellsConnectivityList();
  paraHelper.startCellsoffsetsList();
  
  paraHelper.pushInteger(nb);
  
  paraHelper.endCellsoffsetsList();
  paraHelper.startCellstypesList();

  paraHelper.pushInteger(2);
  
  paraHelper.endCellstypesList();

  if (dumpWeights){
    paraHelper.startPointDataList();
    
    paraHelper.startData("weight",1);
    
    at = itrec.GetFirst();
    unsigned int i = 0;
    while (!itrec.end())
      {
#ifdef DEBUG_BRIDGE
        if (i >= absolute_lambdasA)
          FATAL("depassement detecte");
#endif
        paraHelper.pushDouble(lambdasA[i]);
        at = itrec.GetNext();
        ++i;
      }
    
    paraHelper.endData();
    paraHelper.endPointDataList();
  }

  paraHelper.write_conclusion();

  file.close();

  delete &itrec;

}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::Dump(char * fichier,std::vector<double> & atomes,
                                                    std::vector<int> & suppfield){
  if(suppfield.size()!=atomes.size()/3)
    {
      FATAL("pour la methode dump le champ supplementaire doit etre de la meme taille modulo la dim 3");
    }

  Dump(fichier,atomes,&suppfield[0]);
}

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::Dump(char * fichier,std::vector<double> & atomes,
                                                    int * suppfield){

  UnitsConverter unit;
  unit.SetReadUnits(UnitsConverter::AtomsUnits);

  /***************************************************/
  /* je fait le dump des poids des atomes */
  /***************************************************/

  //  IteratorAtomsRec & itrec = atomes_fictifs.GetIterator();
  if (atomes.size() == 0)
    return;

//   std::vector<double>::iterator itrec = atomes.begin();
//   std::vector<double>::iterator last = atomes.end();

  LMFile file;
  unsigned int nb = atomes.size()/3;
  //unsigned int nb = atomes_fictifs.nbElem();
  
  char filename[256];
  sprintf(filename,"%s-atoms_weigth.vtu",fichier);
  file.open(filename);    

  ParaviewHelper paraHelper;
  paraHelper.setOutputFile(file);
  paraHelper.SetMode(TEXT);

  paraHelper.write_header(nb,1);
  paraHelper.startDofList();
 
  for (unsigned int it=0 ; it < nb ; ++it) 
    {
      double x = atomes[it*3]*unit.etalon_metre;
      double y = atomes[it*3+1]*unit.etalon_metre;
      double z = atomes[it*3+2]*unit.etalon_metre;

      paraHelper.pushDouble(x);
      paraHelper.pushDouble(y);
      paraHelper.pushDouble(z);
    }

  paraHelper.endDofList();
  paraHelper.startCellsConnectivityList();
  
  for (unsigned int i=0;i < nb;++i)
    {
      paraHelper.pushInteger(i);
    }    
  
  paraHelper.endCellsConnectivityList();
  paraHelper.startCellsoffsetsList();
  
  paraHelper.pushInteger(nb);
  
  paraHelper.endCellsoffsetsList();
  paraHelper.startCellstypesList();

  paraHelper.pushInteger(2);
  
  paraHelper.endCellstypesList();

  paraHelper.startPointDataList();

  if (suppfield){
    
    paraHelper.startData("suppfield",1);
    
    for (unsigned int it=0 ; it < nb ; ++it) 
      {
        paraHelper.pushDouble(suppfield[it]);
      }
    
    paraHelper.endData();
  }

  paraHelper.endPointDataList();

  paraHelper.write_conclusion();

  file.close();
}


template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::Dump(char * fichier,
                                                    ContenerNodesRec & noeuds,
                                                    ContenerElemsRec & elems,bool dumpWeights){

  if (elems.nbElem() == 0){
    DUMP("Warning : not dumping to file Bridged information as elems contener seems to be empty",DBG_WARNING);
    return;
  }

  if (noeuds.nbElem() == 0){
    DUMP("Warning : not dumping to file Bridged information as nodes contener seems to be empty",DBG_WARNING);
    return;
  }

  char filename[256];
  sprintf(filename,"%s-mesh_weigth.vtu",fichier);

  LMFile file;

  file.open(filename);    
  
  ParaviewHelper paraHelper;
  paraHelper.setOutputFile(file);
  paraHelper.SetMode(TEXT);

  unsigned int nb=0;

  UnitsConverter unit;
  unit.SetReadUnits(UnitsConverter::AtomsUnits);

  /***************************************************/
  /* je fait le dump des poids du maillage */
  /***************************************************/

  /* definition des types */

  typedef typename DomainC::_Vec_ _Vec_;
  typedef typename DomainC::_Model_ _Model_;
  typedef typename DomainC::RefElt RefElt;
  typedef typename DomainC::ContenerElems ContenerElems;
  typedef typename DomainC::IteratorElems IteratorElems;

  /* structures pour manipuler les noeuds dans le grands vecteur */

  _Model_ & model = static_cast<_Model_ &>(domC.getModel());
  
  _Vec_ & p0 = domC.getDofs().P0();
  nb = noeuds.nbElem();

  /* structures pour manipuler les elements */

  IteratorElemsRec & itc = elems.GetIterator();
  RefElt el = itc.GetFirst();

  std::vector<unsigned int> nodes;
  int compteur;

  /* start */

  paraHelper.write_header(nb,elems.nbElem());
  paraHelper.startDofList();

  for (unsigned int i=0; i < nb; ++i)
    {
      unsigned int n = noeuds.Get(i);

      paraHelper.pushDouble(p0[n*Dim]*unit.etalon_metre);
      if (Dim > 1)
        paraHelper.pushDouble(p0[n*Dim+1]*unit.etalon_metre);
      else
        paraHelper.pushDouble(0.);
      if (Dim > 2)
        paraHelper.pushDouble(p0[n*Dim+2]*unit.etalon_metre);   
      else
        paraHelper.pushDouble(0.);
      
      file.printf("\n");
    }
  
  paraHelper.endDofList();  
  paraHelper.startCellsConnectivityList();

  for (el = itc.GetFirst(); !itc.end() ; el = itc.GetNext())
    {
      model.GlobalIndexes(el,nodes);
      for (unsigned int i=0;i < nodes.size();++i)
        {
          paraHelper.pushInteger(noeuds.Recherche(nodes[i]));
        }
      file.printf("\n");
    }
        
  paraHelper.endCellsConnectivityList();
  
  paraHelper.startCellsoffsetsList();

  compteur = 0;
  for (el = itc.GetFirst(); !itc.end() ; el = itc.GetNext())
    {
      model.GlobalIndexes(el, nodes);
      
      compteur += nodes.size();
      paraHelper.pushInteger(compteur);

      file.printf("\n");
    }

  paraHelper.endCellsoffsetsList();
  paraHelper.startCellstypesList();

  int code_type=-1;
  
  for (el = itc.GetFirst(); !itc.end() ; el = itc.GetNext())
    {
      model.GlobalIndexes(el, nodes);
      switch(nodes.size())
          {
          case 2 : code_type = 3;break;
          case 3 : code_type = 5;break;
          case 8 : code_type = 12;break;
          case 4 : code_type = 10;break;
          default: FATAL("Fatal : unable to dump for VTK : "  << nodes.size());
          }

      paraHelper.pushInteger(code_type);
      file.printf("\n");
    }
  
  paraHelper.endCellstypesList();

  if (dumpWeights){
    paraHelper.startPointDataList();
    paraHelper.startData("weight",1);
    
    for (unsigned int i = 0;i < nb; ++i) 
      {
#ifdef DEBUG_BRIDGE
        if (i >= absolute_lambdasC) FATAL("depassement detecte");
#endif
        paraHelper.pushDouble(lambdasC[i]);
      }
    
    paraHelper.endData();
    paraHelper.endPointDataList();
  }
  
  paraHelper.write_conclusion();
  
  
  file.close();  

  delete &itc;
}

static int cpt=0;

template <typename DomainA, typename DomainC,unsigned int Dim,int Pond>
void Bridging<DomainA,DomainC,Dim,Pond>::DumpAtome(double x,double y, double z){
  //un atome tout seul dans un fichier vtu

  UnitsConverter unit;
  unit.SetReadUnits(UnitsConverter::AtomsUnits);
  

  char filename[255];
  sprintf(filename,"atome-probleme-%d.vtu",cpt);
  ++cpt;

  LMFile file;

  file.open(filename);    

  ParaviewHelper paraHelper;
  paraHelper.setOutputFile(file);
  paraHelper.write_header(1,1);
  paraHelper.startDofList();

  paraHelper.pushDouble(x*unit.etalon_metre);
  paraHelper.pushDouble(y*unit.etalon_metre);
  paraHelper.pushDouble(z*unit.etalon_metre);
  
  paraHelper.endDofList();
  paraHelper.startCellsConnectivityList();

  paraHelper.pushInteger(0);
  
  paraHelper.endCellsConnectivityList();
  paraHelper.startCellsoffsetsList();
  
  paraHelper.pushInteger(1);
  
  paraHelper.endCellsoffsetsList();
  paraHelper.startCellstypesList();

  paraHelper.pushInteger(2);
  
  paraHelper.endCellstypesList();

  paraHelper.write_conclusion();

  file.close();
}


#ifdef USING_STAMP
template class Bridging<DomainStamp,DomainLibMesh<1>,1,LINEAR>;
template class Bridging<DomainStamp,DomainLibMesh<2>,2,LINEAR>;
template class Bridging<DomainStamp,DomainLibMesh<3>,3,LINEAR>;
template class Bridging<DomainStamp,DomainLibMesh<1>,1,CUSTOM>;
#endif 
#ifdef USING_LAMMPS
template class Bridging<DomainLammps,DomainLibMesh<2>,2,LINEAR>;
template class Bridging<DomainLammps,DomainLibMesh<3>,3,LINEAR>;
#endif
template class Bridging<DomainMD1D,DomainLibMesh<1>,1,LINEAR>;
template class Bridging<DomainMD1D,DomainLibMesh<1>,1,CUSTOM>;

---