shape_matrix.h

Go to the documentation of this file.

/* ./bridging/shape_matrix.h 
**********************************
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.
***********************************/

#ifndef SMATRIX_H
#define SMATRIX_H

#include "../math/matrix.h"

template <typename _Vec_,unsigned int Dim>
class ShapeMatrix{
 public:

  ShapeMatrix(int t){
    SubMatrix = new FullMatrix *[t]; 
    tailles = new unsigned int[2*t];
    indirectionAtome[0] = new int *[t];
    indirectionAtome[1] = new int *[t];
    currentIndirectionAtome = indirectionAtome[0];
    newIndirectionAtome = indirectionAtome[1];
    indirectionNode = new int *[t];
    indirectionGlobalNode = new int *[t];
    remplissageA = new int[t];
    remplissageN = new int[t];
    memset(remplissageA,0,sizeof(int)*t);
    memset(remplissageN,0,sizeof(int)*t);

    tab_contigu = NULL;
    taille = t;
  };

  ~ShapeMatrix(){
    for (int i=0;i<taille;++i) 
      delete SubMatrix[i]; 
    
    if (tab_contigu) 
      delete [] tab_contigu; 
    if (indA_contigu[0]) 
      delete [] indA_contigu[0]; 
    if (indA_contigu[1]) 
      delete [] indA_contigu[1]; 
    if (indN_contigu) 
      delete [] indN_contigu; 
    if (indGN_contigu) 
      delete [] indGN_contigu; 
    
    delete [] remplissageA; 
    delete [] remplissageN; 
    delete [] indirectionNode; 
    delete [] indirectionGlobalNode; 
    delete [] indirectionAtome[0]; 
    delete [] indirectionAtome[1]; 
    delete [] SubMatrix; 
    delete [] tailles;
  }

  void SetSub(int i,int m,int n){
    if (m <= 0 || n <= 0)
      FATAL("can not set a smatrix bloc size of zero : m=" << m << " n=" <<n);
    tailles[i*2] = m;
    tailles[i*2+1] = n;
  };

  void Allocate(){
    //premier passage pour calculer tout les
    //coeficients de la smatrice globale
    
    int m,n;

    int total_size=0;
    int total_atoms=0;
    int total_nodes=0;

    for (int i=0;i<taille;++i){
      m = tailles[i*2];
      n = tailles[i*2+1];
      
      total_size+= (m+1)*n;
      total_atoms+=m;
      total_nodes+=n;
    }

    //allocation
    tab_contigu = new double[total_size];
    int index = 0;

    indA_contigu[0] = new int[total_atoms];
    indA_contigu[1] = new int[total_atoms];
    int index_atoms = 0;

    indN_contigu = new int[total_nodes];
    indGN_contigu = new int[total_nodes];
    int index_nodes = 0;

    //deuxieme passage pour allouer les structures de control des matrices
    for (int i=0;i<taille;++i){
      m = tailles[i*2];
      n = tailles[i*2+1];

      SubMatrix[i] = new FullMatrix(m+1,n,tab_contigu+index);

      indirectionAtome[0][i] = indA_contigu[0]+index_atoms;
      indirectionAtome[1][i] = indA_contigu[1]+index_atoms;
      indirectionNode[i] = indN_contigu+index_nodes;
      indirectionGlobalNode[i] = indGN_contigu+index_nodes;

      //incrementation des indexs
      index += (m+1)*n;
      index_atoms += m;
      index_nodes += n;
    }
  }

  void print(){
    for (int k=0;k<taille;++k){
      if (!remplissageA[k]) return;
      FullMatrix & mat = *SubMatrix[k];

      DUMP("affiche les infos de " << k,DBG_DETAIL);

      for (int i =0; i < mat.m() -1; ++i)
        std::cout << indirectionAtome[k][i] << "\t";

      std::cout << endl;
      
      SubMatrix[k]->print();
    }
  }

  void changeAtomIndirection(int el,int old_ind_atome,int new_ind_atome){
    FullMatrix & mat = *SubMatrix[el];
    int  * indA = currentIndirectionAtome[el];
    int  * indAN = newIndirectionAtome[el];
    for (int i = 0 ; i < mat.m()-1;++i)
      if (indA[i] == old_ind_atome) {
        indAN[i] = new_ind_atome;
        break;
      }
  }

  void operator ()(int atome,int node_global,int node,int el,double value){
    int  * indA = currentIndirectionAtome[el];
    indA[remplissageA[el]]=atome;

    int  * indN = indirectionNode[el];
    indN[remplissageN[el]]=node;

    DUMP("setting shapes for atome " << atome << " and element " << el << " node " << node << " node_global " << node_global << " value " << value,DBG_ALL);
    int  * indGN = indirectionGlobalNode[el];
    indGN[remplissageN[el]]=node_global;
    
    // je prend la matrice dense qui correspond a mon element
    FullMatrix & m = *SubMatrix[el];
    m(remplissageA[el],remplissageN[el]) = value;
    // fait les sommes partielles sur les noeuds
    m(m.m()-1,remplissageN[el])+= value;

    ++remplissageN[el];
    if (remplissageN[el] == m.n()){
      remplissageN[el]=0;
      ++remplissageA[el];
    }
  }

  void build_node_shape(double * node_shape){
    for (int k =0 ; k < taille ; ++k){
      int * indN = indirectionNode[k];
      FullMatrix & mat = *SubMatrix[k];
      for (int j = 0 ; j < mat.n(); ++j){
        DUMP("Summing node_shape[" << indN[j] << " (el= " << k << ", j = " << j << ")] with " << mat(mat.m()-1,j),DBG_ALL);
        node_shape[indN[j]] += mat(mat.m()-1,j);
      }
    }
  }

  void buildContinuConstraint(double * A,double * node_shape){
    for (int k =0 ; k < taille ; ++k){
      int * indA = currentIndirectionAtome[k];
      int * indN = indirectionNode[k];
      FullMatrix & mat = *SubMatrix[k];
      for (int i = 0 ; i < mat.m()-1; ++i)
        for (int j = 0 ; j < mat.n(); ++j){
          DUMP("Building A[" << indA[i] << "] with " << mat(i,j) << "*" << node_shape[indN[j]] << " = " 
              << mat(i,j)*node_shape[indN[j]],DBG_ALL);
          A[indA[i]] += mat(i,j)*node_shape[indN[j]];
      }
    }
  }

  void build_rhs(double * rhs,_Vec_ & v){

    //    unsigned int nbT = v.size()/Dim;

    // je parcours tout les blocs    
    for (int k =0 ; k < taille ; ++k){
      //je chope la sous matrice
      FullMatrix & mat = *SubMatrix[k];
      //je chope les indirection sur les noeuds globaux comme 
      // j'ai un acces direct au grand vecteur
      int * indGN = indirectionGlobalNode[k];
      int * indA = currentIndirectionAtome[k];
      // je parcours tout les atomes du bloc
      for (int i = 0 ; i < mat.m()-1; ++i)
        //je parcours tout les noeuds
        for (int j = 0 ; j < mat.n(); ++j){

          DUMP("Before rhs[" << indA[i] << "]=" << rhs[Dim*indA[i]+1],DBG_ALL);

          rhs[Dim*indA[i]] += v[indGN[j]*Dim] * mat(i,j);

          if (Dim > 1)
            rhs[Dim*indA[i]+1] += v[indGN[j]*Dim+1] * mat(i,j);
          if (Dim ==3)
            rhs[Dim*indA[i]+2] += v[indGN[j]*Dim+2] * mat(i,j);

          DUMP("build rhs to " << Dim*indA[i] << " " << Dim*indA[i]+ 1 << " " << Dim*indA[i]+3,DBG_ALL);
          DUMP("Building " << rhs[Dim*indA[i]+1] << "=rhs[" << indA[i] << "] with " << mat(i,j) << "*" << v[indGN[j]*Dim +1] << " = " 
              << v[indGN[j]*Dim+1] * mat(i,j) << " (index = " << indGN[j] << ")",DBG_ALL);

        }
    }
  }


  void apply_changing(_Vec_ & v,double * multL,double * lambdas){
    //    unsigned int nbT = v.size()/Dim;

    // je parcours tout les blocs    
    for (int k =0 ; k < taille ; ++k){
      //je chope la sous matrice
      FullMatrix & mat = *SubMatrix[k];
      //je chope les indirection sur les noeuds globaux comme 
      // j'ai un acces direct au grand vecteur
      int * indGN = indirectionGlobalNode[k];
      int * indN = indirectionNode[k];
      int * indA = currentIndirectionAtome[k];
      // je parcours tout les atomes du bloc
      for (int i = 0 ; i < mat.m()-1; ++i)
        //je parcours tout les noeuds
        for (int j = 0 ; j < mat.n(); ++j){
          DUMP("adding v[" << indGN[j] << "] with " << -1.0*mat(i,j) << " * " <<  multL[indA[i]*Dim+1] << " {rhs[" << indA[i] << "]}/ " 
              << lambdas[indN[j]]<< " = " << -1.0*mat(i,j)*multL[indA[i]*Dim+1]/lambdas[indN[j]],DBG_ALL);

          v.add(indGN[j]*Dim,-1.0*mat(i,j)*multL[indA[i]*Dim]/lambdas[indN[j]]);
          if (Dim > 1)
            v.add(indGN[j]*Dim+1,-1.0*mat(i,j)*multL[indA[i]*Dim+1]/lambdas[indN[j]]);
          if (Dim == 3)
            v.add(indGN[j]*Dim+2,-1.0*mat(i,j)*multL[indA[i]*Dim+2]/lambdas[indN[j]]);
        }
    }
  };

  double InterpolateAtom(int at, _Vec_ & field,int elem,int deg){
    double res = 0;
    int * indA = currentIndirectionAtome[elem];
    int * indGN = indirectionGlobalNode[elem];
    FullMatrix & mat = *SubMatrix[elem];

    DUMP("interpolating for elem " << elem << " atom " << at,DBG_ALL);
    /*     trouve l'atome concerné */
    int i;
    for (i = 0 ; i < mat.m()-1; ++i)
      if (indA[i] == at) 
        break;
    
    /* effectue l'interpolation */

    //    unsigned int nbT = field.size()/Dim;

    for (int j = 0 ; j < mat.n(); ++j){
      res += mat(i,j)*field[indGN[j]*Dim+deg];
    }
    
    return res;
    
  };


  void GetRelatedAtomsToElem(int elem,std::vector<unsigned int> & atoms){
    int * indA = currentIndirectionAtome[elem];
    FullMatrix & mat = *SubMatrix[elem];

    atoms.clear();

    /* parcours les atomes */
    int i;
    for (i = 0 ; i < mat.m()-1; ++i)
      atoms.push_back(indA[i]);
  };


  void GetRelatedNodesToElem(int elem,std::vector<unsigned int> & nodes){
    int * indGN = indirectionGlobalNode[elem];
    FullMatrix & mat = *SubMatrix[elem];

    nodes.clear();

    /* parcours les noeuds */
    int j;
    for (j = 0 ; j < mat.n(); ++j)
      nodes.push_back(indGN[j]);
  };


  /********************* ici c pour least square *****************/

/*   void buildLeastSquareMatrix(FullMatrix * A){ */
  void buildLeastSquareMatrix(FullMatrix * A){
    A->zero();

    /* on veut construire les A(h,k) */

    /* le produit de deux shape sur la position d'un atome
       n'est nul que si les deux shapes font partie
       d'un element commun : on a une seul boucle pricipale sur
       les elements */


    for (int k =0 ; k < taille ; ++k){
      int * indN = indirectionNode[k];
      FullMatrix & mat = *SubMatrix[k];    

      /* ici on met la double boucle sur les noeuds pour calculer la sous matrice de la grande */

      
      for (int j1 = 0 ; j1 < mat.n() ; ++j1)
        for (int j2 = 0 ; j2 < mat.n() ; ++j2)
          {
            /* ici on met la boucle de la somme sur les atomes */
            for (int i = 0 ; i < mat.m()-1 ; ++i){
              (*A)(indN[j1],indN[j2]) += mat(i,j1) * mat(i,j2);
              DUMP("A[" << indN[j1] << "," << indN[j2] << "] = " << (*A)(indN[j1],indN[j2]),DBG_ALL);
            }
          }
    }

    FILE * f = fopen("least-square-constraint.mat","wb+");
    A->PrintToFile(f);
    fclose(f);
  }

  void build_least_square_rhs(Vector * rhs,int at,int elem,double v){
    int * indA = currentIndirectionAtome[elem];
    int * indN = indirectionNode[elem];
    FullMatrix & mat = *SubMatrix[elem]; 
    
    /*     trouve l'atome concerné */
    int i,j;
    for (i = 0 ; i < mat.m()-1; ++i)
      if (indA[i] == at) 
        break;
      
    /* construit le second membre de la methode least square */

    for (j = 0 ; j < mat.n(); ++j)
      (*rhs)(indN[j]) += mat(i,j)*v;

  };
  

  /**************** ici c pour la forme faible **************************/

  void buildFormeFaibleCMatrix(FullMatrix & C,double * lambdas){


    /* le produit de deux shape sur la position d'un atome
       n'est nul que si les deux shapes font partie
       d'un element commun : on a une seul boucle pricipale sur
       les elements */


    for (int k =0 ; k < taille ; ++k){
      int * indN = indirectionNode[k];
      int * indA = currentIndirectionAtome[k];

      FullMatrix & mat = *SubMatrix[k];    

      /* ici on met la double boucle sur les noeuds pour calculer la sous matrice de la grande */

      for (int j1 = 0 ; j1 < mat.n() ; ++j1)
        for (int j2 = 0 ; j2 < mat.n() ; ++j2)
          {
            /* ici on met la boucle de la somme sur les atomes */
            for (int i = 0 ; i < mat.m()-1 ; ++i){

              C(indN[j1],indN[j2]) += mat(i,j1) * mat(i,j2) / lambdas[indA[i]];
              DUMP("C[" << indN[j1] << "," << indN[j2] << "] = " << C(indN[j1],indN[j2]),DBG_ALL);
            }
          }
    }
    
  };



  void correct_forme_faible(double & v,Vector & rhs,double * lamb,int at,unsigned int elem){
    int * indN = indirectionNode[elem];
    int * indA = currentIndirectionAtome[elem];
    FullMatrix & mat = *SubMatrix[elem];    

    /*     trouve l'atome concerné */
    int i,j;
    for (i = 0 ; i < mat.m()-1; ++i)
      if (indA[i] == at) 
        break;


    for (j = 0 ; j < mat.n() ; ++j){
      v -= mat(i,j)*rhs(indN[j])/lamb[at];
    }

  };

  //****************************************************************
  //network operations
  //****************************************************************

/*   void TransitThroughNetwork(Communicator & communicator,DuoDistributedVecteur & coms,int groupefrom,int groupedest){  */
/*     // premiere chose, je fait une boucle sur les blocks car c'est eux qui conditionnent le nombre d'envois  */
/*     unsigned int nb_coms = coms.nb_blocks();  */

/*      DUMP("on va faire " << nb_coms << " communications");  */
    
/*      for(unsigned int com=0; com < nb_coms ; ++com){  */
/*        int offset=0;  */
/*        int sendsize=0;  */
    
/*        // comme les atomes sont rangés par elements et non par atomes on va devoir faire des parcours pour filtrage  */
/*        //on commence le parcours   */
      
/*        for (int elem = 0 ; elem < taille ; ++elem){    */
/*      int nb_atoms_for_com=0;  */
        
/*      int * indA = indirectionAtome[elem];  */
/*      int m = tailles[elem*2]; */
/*      int n = tailles[elem*2+1]; */

/*      //je compte les atomes qui sont pour com  */
/*      for (int i = 0 ; i < m;++i)  */
/*        if(coms.isInBlock(indA[i],com))  */
/*          ++nb_atoms_for_com;  */

/*      //premier cas : la matrice ne contient que des atomes pour le destinataire com   */
/*      //je continue je n'envoi pas de suite  */
/*      if(nb_atoms_for_com==m){  */
/*        //j'augmente la taille a envoyer de 1 element  */
/*        ++sendsize; */
/*        continue; */
/*      }  */

/*      FullMatrix & mat = *SubMatrix[offset];  */

/*      //dans les autres cas je commence par envoyer ce qui est en attente  */
/*      //pour cela je commence par le tableau de tailles qui correspond aux matrices en attentes  */
/*      communicator.SendIntegers(&sendsize,1,com,groupedest);  */
/*      communicator.SendIntegers(tailles+2*offset,sendsize*2,com,groupedest);  */
/*      //j'envoie ensuite les vecteur d'indirection  */
/*      communicator.SendIntegers(indirectionAtome[offset], */
/*                                indirectionAtome[offset+sendsize+1]-indirectionAtome[offset], */
/*                                com,groupedest);  */
/*      communicator.SendIntegers(indirectionNode[offset], */
/*                                indirectionNode[offset+sendsize+1]-indirectionNode[offset], */
/*                                com,groupedest);  */
/*      //finalement j'envoie le gros bloc des resultats de shapes */
/*      commnuicator.SendDoubles(&mat(0,0)); */

/*      //je remet a zero sendsize et je remet l'offset sur l'element courant */
/*      offset = elem; */
/*      sendsize = 0; */

/*      //deuxieme cas : la matrice ne contient aucun atomes pour le destinataire com */
/*      // je poursuit la boucle */
/*      if(nb_atoms_for_com==0){ */
/*        ++offset; */
/*        continue; */
/*      } */

/*      //troisieme cas : la matrice contient certain atomes pour le destinataire com */
/*      else{ */
/*        sendsize=1; */

/*        //alors je reconstruit une matrice que j'envoie */
/*        int new_m = nb_atoms_for_com+1; */
/*        double * mat_tmp = new double[new_m*n]; */
/*        memset(mat_tmp,0,new_m*n); */

/*        int taille_tmp[2]; */
/*        taille_tmp[0] = mew_m; */
/*        taille_tmp[1] = n; */
          
/*        int index = 0; */

/*        //je rempli la nouvelle matrice */
/*        for (int i = 0 ; i < m;++i)  */
/*          if(coms.isInBlock(indA[i],com)){ */
/*            for (int j=0; j < n ; ++j){ */
/*              mat_tmp[index+new_m*j] = mat(i,j); */
/*              mat_tmp[new_m-1+new_m*j]+=mat(i,j); */
/*            } */
/*            ++index; */
/*          } */
          
/*        //je l'envoie de la meme maniere par le reseau */
/*        communicator.SendIntegers(&sendsize,1,com,groupedest); */
/*        communicator.SendIntegers(taille_tmp,2,com,groupedest); */
/*        //j'envoie les vecteurs d'indirection */

/*        delete [] mat_tmp; */
          
/*      } */
          
          
/*        } */
/*      } */
/*   } */
    
  int * indirectionAtomTable(){
    return indA_contigu;
  }

  void swapAtomIndirections(){
    int ** temp = currentIndirectionAtome;
    currentIndirectionAtome = newIndirectionAtome;
    newIndirectionAtome = temp;
  }
  

 private:

  //  friend class Bridging;  
  int taille;

  unsigned int * tailles;

  FullMatrix ** SubMatrix;
  int ** currentIndirectionAtome; // indirection vers le vecteur pricipal atomes
  int ** newIndirectionAtome; // indirection vers le vecteur pricipal atomes

  int ** indirectionAtome[2]; // indirection vers le vecteur pricipal atomes

  int ** indirectionNode; // indirection vers le vecteur principal noeuds
  int ** indirectionGlobalNode; // indirection vers des references de noeuds
  int * remplissageA; // indice de remplissage de la matrice locale selon les atomes
  int * remplissageN; // indice de remplissage de la matrice locale selon les noeuds

  //ici je place les tableau contigu des données pour etre sur que la structure est bien
  //contigu je fait des allocation spéciales
  
  double * tab_contigu;  
  int * indA_contigu[2];

  int * indN_contigu;
  int * indGN_contigu;
};

#endif

---