BoxPartition.cpp

Go to the documentation of this file.
00001 // This is a -*- C++ -*-  file
00002 //------------------------------------------------------------------------------------//
00003 //                                                                                    //
00004 //   Copyright (c) 2001                                                               //
00005 //      INRIA                                                                         //
00006 //      54600 VILLERS LES NANCY                                                       //
00007 //      France                                                                        //
00008 //                                                                                    //
00009 //------------------------------------------------------------------------------------//
00010 //                                                                                    //
00011 //                 * NOTICE OF PROPRIETARY INFORMATION *                              //
00012 //                                                                                    //
00013 // The information contained in this file is considered proprietary and the exclusive //
00014 // property of  INRIA. This information may not be disclosed, duplicated              //
00015 // or used, in whole or in part, for  any purpose  whatsoever without express written //
00016 // authorization from INRIA                                                           //
00017 //     
00018 
00019 #include "BoxPartition.h"
00020 #if defined (_GNU_GSL_)
00021 #include <gsl/gsl_math.h>
00022 #include <gsl/gsl_eigen.h>
00023 #else
00024 #include "newTools.h"
00025 #endif
00026 //
00027 #include <iomanip>
00028 #include <vector>
00029 //
00030 BoxPartition::BoxPartition() : PartitionBase< setOfAtoms >() 
00031 {
00032   _gridOfPartition[0] = _gridOfPartition[1] = _gridOfPartition[2] = 0;
00033   this->setToZero() ; 
00034 }
00035 
00036 BoxPartition::BoxPartition(const int size, const double &R1, const double & R2) : PartitionBase< setOfAtoms >(size,R1,R2)
00037 {
00038   _gridOfPartition[0] = size ; _gridOfPartition[1] = _gridOfPartition[2] = 1;
00039   this->setToZero() ;
00040 }
00041 
00042 BoxPartition::BoxPartition(const int size1, const int size2, const int size3, const double &R1, const double & R2) : 
00043                  PartitionBase< setOfAtoms >(size1*size2*size3,R1,R2)
00044 {
00045   _gridOfPartition[0] = size1 ; _gridOfPartition[1] = size2 ; _gridOfPartition[2] = size3;
00046 }
00047   //
00048 BoxPartition::BoxPartition(const int size[3], const double &R1, const double & R2) :  PartitionBase< setOfAtoms >( size[0] * size[1] *  size[2],R1,R2)
00049 {
00050   _gridOfPartition[0] = size[0] ; _gridOfPartition[1] = size[1] ; _gridOfPartition[2] =  size[2] ;
00051 }
00052 
00053 BoxPartition::~BoxPartition()
00054 { this->setToZero() ; }
00055 
00056 void BoxPartition::setToZero()
00057 {
00058   for( int i = 0 ; i < 3 ; ++i) {
00059     _gravityCenter[i]     = 0.0 ; 
00060     _axe1[i]              = 0.0 ;
00061     _axe2[i]              = 0.0 ;
00062     _axe3[i]              = 0.0 ;
00063     _inertialMatrix[0][i] = 0.0 ;
00064     _inertialMatrix[1][i] = 0.0 ;
00065     _inertialMatrix[2][i] = 0.0 ;
00066   }
00067 }
00068 void BoxPartition::innerReperOfMolecule(){
00069 
00070   std::cout << "  BoxPartition::innerReperOfMolecule : " <<std::endl
00071             << "      nbAtoms "<< _setOfAtoms->getNumberOfAtoms() <<std::endl;
00072   _setOfAtoms->buildGravityCenter(_gravityCenter) ;
00073   std::cout<<"Centre de Gravite : ("<<_gravityCenter[0]<<" ; "<<_gravityCenter[1]<<" ; "
00074            <<_gravityCenter[2]<<")"<<std::endl;
00075   //
00076   int i ;
00077   double Ixx, Iyy, Izz, Ixy, Ixz, Iyz ;
00078   Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0 ;
00079   for(i = 0 ; i < _setOfAtoms->getNumberOfAtoms() ; ++i)
00080     {
00081       atomElement & atoms = this->_setOfAtoms->getElement(i) ;
00082       //      std::cout << i  << " pos " <<  atoms._name<< " "<<   atoms._xyz[1] <<std::endl ;
00083       Ixx += ( atoms._xyz[1] - _gravityCenter[1])*( atoms._xyz[1] - _gravityCenter[1]) +
00084              ( atoms._xyz[2] - _gravityCenter[2])*( atoms._xyz[2] - _gravityCenter[2]);
00085       Iyy += ( atoms._xyz[0] - _gravityCenter[0])*( atoms._xyz[0] - _gravityCenter[0]) +
00086              ( atoms._xyz[2] - _gravityCenter[2])*( atoms._xyz[2] - _gravityCenter[2]);
00087       Izz += ( atoms._xyz[0] - _gravityCenter[0])*( atoms._xyz[0] - _gravityCenter[0]) +
00088              ( atoms._xyz[1] - _gravityCenter[1])*( atoms._xyz[1] - _gravityCenter[1]);
00089       Ixy -= ( atoms._xyz[0] - _gravityCenter[0])*( atoms._xyz[1] - _gravityCenter[1]);
00090       Ixz -= ( atoms._xyz[0] - _gravityCenter[0])*( atoms._xyz[2] - _gravityCenter[2]);
00091       Iyz -= ( atoms._xyz[2] - _gravityCenter[2])*( atoms._xyz[1] - _gravityCenter[1]);
00092     }//for i
00093   //
00094   _inertialMatrix[0][0] = Ixx ;
00095   _inertialMatrix[1][1] = Iyy ;
00096   _inertialMatrix[2][2] = Izz ;
00097   _inertialMatrix[0][1] = _inertialMatrix[1][0] = Ixy ;
00098   _inertialMatrix[0][2] = _inertialMatrix[2][0] = Ixz ;
00099   _inertialMatrix[1][2] = _inertialMatrix[2][1] = Iyz ;
00100   double eigenValues[3] ;
00101   int nbRoots = 0 ;
00102 
00103 #if defined (_GNU_GSL_)
00104   gsl_matrix_view m 
00105     = gsl_matrix_view_array(&_inertialMatrix[0][0], 3, 3);
00106   gsl_vector *eval = gsl_vector_alloc (3);
00107   gsl_matrix *evec = gsl_matrix_alloc (3, 3);
00108 
00109   gsl_eigen_symmv_workspace * w = 
00110     gsl_eigen_symmv_alloc (3);
00111   
00112   gsl_eigen_symmv (&m.matrix, eval, evec, w);
00113 
00114   gsl_eigen_symmv_free(w);
00115 
00116   gsl_eigen_symmv_sort (eval, evec, 
00117                         GSL_EIGEN_SORT_VAL_ASC);
00118    for(i = 0 ; i < 3 ; ++i ){
00119      eigenValues[i] = gsl_vector_get(eval, i) ;
00120    }
00121    int j ;
00122   gsl_vector_view evec_i = gsl_matrix_column(evec, 0);
00123   for ( j = 0 ; j < 3 ; ++j){
00124     _axe1[j] = gsl_vector_get(&evec_i.vector, j) ;
00125     //    _axe1[j] = gsl_matrix_get(m, i, j); gsl_vector_get(evec_i, j) ;
00126   }
00127    evec_i = gsl_matrix_column(evec, 1);
00128    for ( j = 0 ; j < 3 ; ++j){
00129      _axe2[j] = gsl_vector_get(&evec_i.vector, j) ;
00130    }
00131    evec_i = gsl_matrix_column(evec, 2);
00132    for( j = 0 ; j < 3 ; ++j){
00133      _axe3[j] = gsl_vector_get(&evec_i.vector, j) ;
00134    }
00135 #else
00136    //
00137    //#error("rien a faire ici")
00138 //   std::cout << std::setprecision(16) << setiosflags(ios::scientific);
00139 //   for(i = 0 ; i < 3 ;++i){
00140 //     for(int j = 0 ; j < 3 ;++j){
00141 //       std::cout << "  "<< _inertialMatrix[i][j] ;
00142 //     }
00143 //     std::cout << std::endl;
00144 //   }
00145 //   std::cout << std::setprecision(6);
00146 //   double coeffs[] =  { -1, Ixx+Iyy+Izz,(-Ixx*Iyy-Iyy*Izz-Ixx*Izz+Ixy*Ixy+Iyz*Iyz+Ixz*Ixz),
00147 //                     Ixx*Iyy*Izz-Ixx*Iyz*Iyz-Ixy*Ixy*Izz+2*Ixy*Ixz*Iyz-Ixz*Ixz*Iyy,
00148 //                       };
00149 //   std::cout << "Coeficients : "<< std::endl
00150 //          << coeffs[0] << "  "<<coeffs[1] << "  "<<coeffs[2] << "  "<<coeffs[3] << "  "<< std::endl;
00151   eigenValues[0] = eigenValues[1] = eigenValues[2] = 0.0 ;
00152   findRootsOfPolynomOfDegre3(coeffs, nbRoots ,eigenValues) ;
00153   double val ; 
00154   for(i = 0 ; i < 3 ; ++i ){
00155     val =coeffs[3]+ eigenValues[i]*( coeffs[2]+ eigenValues[i]*( coeffs[1]+ eigenValues[i]* coeffs[0])) ;
00156     //    std::cout << i << "  "<< eigenValues[i] << "  "<< val <<std::endl;
00157   }
00158   //
00159   //  std::cout << "Calcul des Vecteurs propres " <<std::endl;
00160   findEigenVector(_inertialMatrix,eigenValues[0], _axe1) ;
00161   findEigenVector(_inertialMatrix,eigenValues[1], _axe2) ;
00162   findEigenVector(_inertialMatrix,eigenValues[2], _axe3) ;
00163 #endif
00164   //
00165   std::cout << "Les valeurs Propres et Vecteurs propres : "<< nbRoots<<std::endl
00166             << "    "<< eigenValues[0]<< "  "<<_axe1[0]<< "  "<< _axe1[1]<< "  "<< _axe1[2]<<std::endl
00167             << "    "<< eigenValues[1]<< "  "<<_axe2[0]<< "  "<< _axe2[1]<< "  "<< _axe2[2]<<std::endl
00168             << "    "<< eigenValues[2]<< "  "<<_axe3[0]<< "  "<< _axe3[1]<< "  "<< _axe3[2]<<std::endl
00169             << std::endl;
00170     //
00171 
00172 }
00173 void BoxPartition::projectOnInnerReper(){
00174   //
00175   double s1,s2,s3 ;
00176   s1 =  this->_setOfAtoms->getElement(0)._xyz[0]*_axe1[0] +  
00177     this->_setOfAtoms->getElement(0)._xyz[1]*_axe1[1] +  
00178     this->_setOfAtoms->getElement(0)._xyz[2]*_axe1[2] ;
00179 //   if(s1 > 0) 
00180 //     { std::cout << "Il y a un probleme le repere est mal oriente" << std::endl;
00181 //     _axe1[0] = -_axe1[0] ;_axe1[1] = -_axe1[1] ;_axe1[2] = -_axe1[2] ;
00182 //     }
00183   for(int i = 0 ; i < _setOfAtoms->getNumberOfAtoms() ; ++i)
00184     {
00185       atomElement & atoms = this->_setOfAtoms->getElement(i) ;
00186       atoms._xyz[0] -= _gravityCenter[0];
00187       atoms._xyz[1] -= _gravityCenter[1];
00188       atoms._xyz[2] -= _gravityCenter[2];
00189       s1 =  atoms._xyz[0]*_axe1[0] +  atoms._xyz[1]*_axe1[1] +  atoms._xyz[2]*_axe1[2] ;
00190       s2 =  atoms._xyz[0]*_axe2[0] +  atoms._xyz[1]*_axe2[1] +  atoms._xyz[2]*_axe2[2] ;
00191       s3 =  atoms._xyz[0]*_axe3[0] +  atoms._xyz[1]*_axe3[1] +  atoms._xyz[2]*_axe3[2] ;
00192       atoms._xyz[0] = s1 ;
00193       atoms._xyz[1] = s2 ;
00194       atoms._xyz[2] = s3 ;
00195       //      std::cout << "  "<<i << "   " << s1 << "  " << s2 << "  " << s3 << "  " << _charge[i] <<std::endl;
00196     }
00197 }
00198 bool  BoxPartition::partitionate(int numberOfIterations,  const bool use_frag){
00199   //
00200   const double &R_overlap = _Rbuffer1 ; 
00201   const double &R_cut     = _Rbuffer2 ;
00202   int numberOfAtoms = this->_setOfAtoms->getNumberOfAtoms() ;
00203   //
00204   std::cout << "partitionate " << numberOfAtoms <<std::endl;
00205   //
00206   int i ;
00207   double LengthOfMolecule[3], coinMin[3], coinMax[3] ;
00208   double sizeOfKernel[3], s1[3] , s2[3] ;
00209   bool isCV ;
00210   _useFrag = use_frag ;
00211   //
00212   //  On se place dans le repere interne de la molecule
00213   //
00214   this->innerReperOfMolecule() ;
00215   this->projectOnInnerReper() ;
00216   //
00217   //  Trier les atomes dans ce nouveau repere
00218   //
00219   int* ordre ; ordre = new int [numberOfAtoms] ;
00220   for(i = 0 ; i < numberOfAtoms ; ++i) {
00221     ordre[i] = i;
00222   }//for i
00223   this->sortInFirstDimension(ordre) ; //version damien
00224   // 
00225   //  Structure des voisins
00226   //
00227   //  std::cout << "  ---> struture des voisins " << std::endl;
00228   this->setNeigboors();
00229   //
00230   if(_useFrag){
00231     build_fragmentStructure( ) ;
00232   }
00233   //
00234   //  Init the size of each subdomain
00235   //
00236   //  std::cout << "  ---> find minmax " << std::endl;
00237   this->_setOfAtoms->findMinMax(coinMin , coinMax  ) ;
00238   //
00239   LengthOfMolecule[0]  =   coinMax[0] - coinMin[0] ;
00240   LengthOfMolecule[1]  =   coinMax[1] - coinMin[1] ;
00241   LengthOfMolecule[2]  =   coinMax[2] - coinMin[2] ;
00242   //
00243   //
00244   std::cout << " XMIN            : "<< coinMin[0] << "  "<< coinMin[1] << "  "<< coinMin[2] << std::endl;
00245   std::cout << " XMAX            : "<< coinMax[0] << "  "<< coinMax[1] << "  "<< coinMax[2] << std::endl;
00246   std::cout << "LengthOfMolecule : " <<LengthOfMolecule[0]<< "  "<<LengthOfMolecule[1]
00247             << "  "<<LengthOfMolecule[2]<< "  "<<std::endl;
00248   //
00249   //     Uniquement le cas lineaire
00250   //
00251   sizeOfKernel[0] = LengthOfMolecule[0] /static_cast<double>(_gridOfPartition[0]) ;
00252   sizeOfKernel[1] = LengthOfMolecule[1] /static_cast<double>(_gridOfPartition[1]); 
00253   sizeOfKernel[2] = LengthOfMolecule[2] /static_cast<double>(_gridOfPartition[2]);
00254   //
00255   //
00256   //
00257   for( i = 0 ; i < 3 ; ++i){
00258     //    std::cout << "_gridOfPartition : "<<_gridOfPartition[i] << " sizeOfKernel : " << sizeOfKernel[i]<<std::endl;
00259 //   s1[i]            = coinMin[i] ;
00260 //     s2[i]            = coinMin[i] + sizeOfKernel[i];
00261   }
00262   //  std::cout << "  ---> Set min and max in each subdomain " << std::endl;
00263   int nyz , nsd ;
00264   double coin1[3], coin2[3];
00265   for( i = 0 ; i < 3 ; ++i){
00266     coin1[i] = coinMin[i] ;
00267     coin2[i] = coinMax[i] ;
00268     s1[i]    = coinMin[i] ;
00269     s2[i]    = coinMax[i] ;
00270   }
00271   for (int nsdz = 0 ; nsdz < _gridOfPartition[2]  ; ++nsdz){
00272     coin1[2]  = coinMin[2] + static_cast<double>( nsdz ) * sizeOfKernel[2];
00273     coin2[2]  = coin1[2] + sizeOfKernel[2];
00274     s1[2]     =  coin1[2] ;
00275     s2[2]     =  coin2[2] ;
00276     for (int nsdy = 0 ; nsdy < _gridOfPartition[1] ; ++nsdy){
00277       nyz = _gridOfPartition[0]*(nsdy+_gridOfPartition[1]* nsdz);
00278       s1[1]  = coin1[1] + static_cast<double>( nsdy  ) * sizeOfKernel[1];
00279       s2[1]  = s1[1] + sizeOfKernel[1];
00280       //
00281       s1[0]  = coinMin[0] ;
00282       s2[0]  = coinMin[0] + static_cast<double>( _gridOfPartition[0]) * sizeOfKernel[0];
00283       _subDomains[nyz].setXmin(s1) ; 
00284       _subDomains[_gridOfPartition[0]-1+ nyz].setXmax(s2) ; 
00285       for (int nsdx = 1 ; nsdx < _gridOfPartition[0] ; ++nsdx){
00286         nsd   = nsdx +  nyz ;
00287         s1[0] = coin1[0] + static_cast<double>( nsdx) * sizeOfKernel[0];
00288         s2[0] =  s1[0] ;//+  sizeOfKernel[0];
00289 //      std::cout <<  nsd << " "<< nsdx << " "<< nsdy << " "<< nsdz << " s1 s2 coin1 coin2 " << std::endl;
00290 //      for (int q = 0 ; q<3 ; ++q){
00291 //        std::cout <<"                          "<< s1[q] << "  "<< s2[q] << "  "<< coin1[q] << "  "<< coin2[q] << "  "<<std::endl;
00292 //      }
00293         _subDomains[nsd-1].setXmax( s2)  ; 
00294         _subDomains[nsd].setXmin( s1); 
00295       }
00296     }
00297   }
00298   for (int nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00299     std::cout << "=====> Length : " <<  _subDomains[nsd].getXLength()  << "  "
00300               <<  _subDomains[nsd].getYLength()  << "  "<<  _subDomains[nsd].getZLength()  << std::endl;
00301   } 
00302 
00303   // ------------------------------------------------------------
00304   // Boucle pour l'équilibrage du nombre d'OA par SD :
00305   // ------------------------------------------------------------
00306   //
00307   int equilibrateCpt;
00308   int nbMaxOfIterationToEquilibrate = numberOfIterations;
00309   for(equilibrateCpt = 0 ; equilibrateCpt < nbMaxOfIterationToEquilibrate ; ++equilibrateCpt) {
00310     //
00311     // Add atoms in the kernel
00312     //
00313     //    std::cout << "  ---> set Atoms In Kernel " << std::endl;
00314     this->setAtomsInKernel() ;
00315     //
00316     //
00317     // Add atoms in Overlap and Cut
00318     //
00319     //    std::cout << "  ---> set Atoms In Overlap and Cut " << R_overlap << "  "<< R_cut <<std::endl;
00320     //   this->setAtomsInOverlapAndCut(R_overlap, R_cur); 
00321     this->setAtomsInOverlapAndCut(); 
00322     //
00323     //    std::cout << "  ---> set numbre Of Atoms in subdomains  " <<std::endl;
00324     for (int nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00325       _subDomains[nsd].computeNumberOfAtoms();
00326     }
00327     //    std::cout << "  ---> set Load of subdomains "  <<std::endl;
00328     this->setLoad() ;
00329     //
00330     //    std::cout << "  ---> computeBoxEquilibrateCharge "  <<std::endl;
00331     std::cout <<" Iteration : " << equilibrateCpt <<std::endl;
00332     isCV = computeBoxEquilibrateCharge(R_overlap,  R_cut);
00333     if(isCV || nbMaxOfIterationToEquilibrate == equilibrateCpt + 1)
00334       {break ;}
00335     // resetData
00336     //    std::cout << "  ---> resetData "  <<std::endl;
00337     this->resetData();
00338   }
00339   //
00340   //  std::cout << "  ---> buildLocalNumerotation "  <<std::endl;
00341   this->buildLocalNumerotation();
00342   //
00343   // For QCPP
00344   //
00345   this->buildAtomsInZone() ;
00346   //
00347   //
00348  return true ;
00349 }
00350 
00351 void BoxPartition::setNeigboors(){
00352   //
00353   int nsd ;
00354   //
00355   int nsdx, nsdy, nsdz, Nyz, Nz ;
00356   //
00357   // Set Global Number
00358   //
00359   //  std::cout << "Set global number OK" << std::endl;
00360   for ( nsdz = 0 ; nsdz <  _gridOfPartition[2] ; ++nsdz){
00361     for ( nsdy = 0 ; nsdy <  _gridOfPartition[1] ; ++nsdy){
00362       for ( nsdx = 0 ; nsdx <  _gridOfPartition[0] ; ++nsdx){
00363         nsd                         = nsdx +  _gridOfPartition[0]*( nsdy +  _gridOfPartition[1]*nsdz );
00364         //      _subDomains[nsd].myNumber() = nsd ;
00365         _subDomains[nsd].setMyNumber(nsd,nsdx ,nsdy, nsdz);
00366       }
00367     }
00368   }
00369   //
00370   // Set Neigboors 
00371   //
00372   int  kz, kz_deb, kz_fin, ky, ky_deb, ky_fin, kx, kx_deb, kx_fin ;
00373   //  std::cout << "Set in x " << std::endl;
00374   for ( nsdz = 0 ; nsdz < _gridOfPartition[2] ; ++nsdz){
00375     if( _gridOfPartition[2]!= 1){
00376       if( nsdz ==0)
00377         { kz_deb = 0; kz_fin = 1; }
00378       else if( nsdz == _gridOfPartition[2] -1)
00379         { kz_deb = -1; kz_fin = 0; }
00380       else {
00381         kz_deb = -1; kz_fin = 1; 
00382       }
00383     }
00384     else {
00385       kz_deb = kz_fin = 0 ;
00386     }
00387     Nz = _gridOfPartition[1]*nsdz ;
00388     for( kz = kz_deb ; kz <= kz_fin ; ++kz){
00389       Nz = _gridOfPartition[1]*( nsdz + kz ) ;
00390       for ( nsdy = 0 ; nsdy < _gridOfPartition[1] ; ++nsdy){
00391         if( _gridOfPartition[1]!= 1){
00392           if( nsdy ==0)
00393             { ky_deb = 0; ky_fin = 1; }
00394           else if( nsdy == _gridOfPartition[1] -1)
00395             { ky_deb = -1; ky_fin = 0; }
00396           else {
00397             ky_deb = -1; ky_fin = 1; }
00398         }       
00399         else {
00400           ky_deb = ky_fin = 0 ;
00401         }
00402         for( ky = ky_deb ; ky <= ky_fin ; ++ky){
00403           Nyz = _gridOfPartition[0]*( nsdy +ky +  Nz) ;
00404           
00405           for ( nsdx = 0 ; nsdx < _gridOfPartition[0] ; ++nsdx){
00406             nsd = nsdx +  _gridOfPartition[0]*( nsdy +_gridOfPartition[1]* nsdz) ;
00407               if( _gridOfPartition[0]!= 1){
00408                 if( nsdx ==0){ 
00409                   kx_deb = 0 ; kx_fin = 1; }
00410                 else if( nsdx == _gridOfPartition[0] -1)
00411                 { kx_deb = -1; kx_fin = 0; }
00412                 else {
00413                   kx_deb = -1; kx_fin = 1; 
00414                 }
00415               } 
00416               else {
00417                 kx_deb = kx_fin = 0 ;
00418               }
00419             for( kx = kx_deb ; kx <= kx_fin ; kx++){
00420               // std::cout << "        nsdx : "<<  nsdx << " "<< nsdy <<" "<<nsdz <<"  kx ky kz " << kx <<" "<<ky<< " "<< kz<<std::endl;
00421               if(kx==0 && ky==0 && kz==0 ){
00422                 continue;}
00423               _subDomains[nsd].addNeigboors(nsdx + kx + Nyz );
00424             }
00425           }
00426         }
00427       }
00428     }
00429   }
00430   //
00431   //  Allocation
00432   //
00433   for ( nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00434     _subDomains[nsd].initNumberOfOverlapDomain() ;
00435   }
00436   //
00437   // Affichage 
00438   //
00439   for ( int nsdz = 0 ; nsdz <  _gridOfPartition[2] ; ++nsdz){
00440     for ( int nsdy = 0 ; nsdy <  _gridOfPartition[1] ; ++nsdy){
00441       for ( int  nsdx = 0 ; nsdx <  _gridOfPartition[0] ; ++nsdx){
00442         nsd = nsdx +  _gridOfPartition[0]*( nsdy +  _gridOfPartition[1]*nsdz );
00443         //      std::cout << " Gnum : "<<_subDomains[nsd].myNumber() << " - "<< nsdx<< " "<< nsdy<< " "
00444         //                << nsdz <<  "  "<< _subDomains[nsd].getNumberOfNeighbors() << std::endl 
00445 //                << "            " ;
00446 //      for (int k = 0; k < _subDomains[nsd].getNumberOfNeighbors() ; ++k){
00447 //        std::cout << "  "<<  _subDomains[nsd].getNumberOfNeighbor(k) ;
00448 //      }
00449 //      std::cout << std::endl;
00450       }
00451     }
00452   } 
00453   //
00454 }
00455 bool  BoxPartition::setAtomsInKernel(){
00456 
00457  bool retVal = false;
00458  int numberOfAtoms = this->_setOfAtoms->getNumberOfAtoms(), nsd, j ;
00459  std::vector<int> emptyKernel;
00460  //
00461  // std::cout << "BoxPartition::setAtomsInKernel() "  <<_numberOfPartition<<std::endl;
00462  //
00463  j = 0 ;
00464  for ( nsd = 0 ; nsd < _numberOfPartition; ++nsd){
00465    //   _subDomains[nsd].addAtomsInKernel(*_setOfAtoms, _useFrag);
00466    //   void subDomain::addAtomsInKernel( const setOfAtoms &atoms, bool use_frag){
00467    //
00468    // Uniquement valable pour le 1D
00469      double size[3] ,dist[3], xmin[3];
00470      AtomIn  atomNSD ;
00471      //
00472      _subDomains[nsd].getLength(&(size[0]));
00473      _subDomains[nsd].getXmin(&(xmin[0]));
00474      //   std::cout << _myNum << "  xmin : "<<  _xmin[0]<<  "  "<<  _xmin[1]<<  "  "<<  _xmin[2]
00475      //             << "  xmax : "<<  _xmax[0]<<  "  "<<  _xmax[1]<<  "  "<<  _xmax[2] << std::endl;
00476      atomNSD._numDomain = nsd;  atomNSD._typeZone = Core_Zone ; atomNSD._localNum = -1 ;
00477      int j = 0 ;
00478      //
00479      for (j = 0 ; j < _setOfAtoms->getNumberOfAtoms() ; ++j){
00480        //
00481        dist[0]     = _setOfAtoms->getElement(j)._xyz[0] - xmin[0] ;
00482        dist[1]     = _setOfAtoms->getElement(j)._xyz[1] - xmin[1] ;
00483        dist[2]     = _setOfAtoms->getElement(j)._xyz[2] - xmin[2] ;
00484        //
00485        if ( ((dist[0] >= 0 && dist[0] <= size[0]) &&  (dist[1] >= 0 && dist[1] <= size[1]) &&  (dist[2] >= 0 && dist[2] <= size[2]) )) {
00486          //      cout<<"L'atome "<< j <<" est dans le SD "<<_myNum<< "  " <<atoms.getElement(j)._xyz[0]-xmin<<"  "<<size<< endl;
00487          _subDomains[nsd].addIndexInKernel(j) ;
00488        }
00489        //    }
00490      }//dfr
00491      std::vector<int> & atomIsInKernel =  _subDomains[nsd].getAtomsIsInKernel() ;
00492      std::vector<int> & indexInKernel  =  _subDomains[nsd].getAtomsIndexInKernel() ;
00493      atomIsInKernel.resize(indexInKernel.size());
00494      for(j= 0 ; j< static_cast<int>(indexInKernel.size()) ; ++j){
00495        atomIsInKernel[j] = 1;
00496      } 
00497      //
00498  }
00499 
00500  int nbAtomsInSubDomains = 0 , num;
00501  for ( nsd = 0 ; nsd < _numberOfPartition; ++nsd){
00502    num                  =_subDomains[nsd].numberOfAtomsInKernel() ;
00503    nbAtomsInSubDomains += num ;
00504    if (num ==0)
00505      {emptyKernel.push_back(nsd) ;
00506      std::cerr << "Pas d'atomes dans le kernel du  sousdomaine "<<nsd<<std::endl;}
00507    //   std::cout << nsd << "  " << _subDomains[nsd].getXLength() << "  "<< num  << std::endl ;
00508  }
00509  if(numberOfAtoms-nbAtomsInSubDomains!=0)
00510    {std::cerr << "NumberOfAtoms : "<< numberOfAtoms<< "  Perte d'atomes : "<< numberOfAtoms-nbAtomsInSubDomains<<std::endl;}
00511  //
00512  // On enleve les domaines vides
00513  //
00514  if(emptyKernel.size() !=0 )
00515    {
00516      for ( nsd = 0 ; nsd < _numberOfPartition; ++nsd){
00517        //       _subDomains[nsd].removeNeigboors(emptyKernel) ;
00518      }
00519    }
00520  return retVal ;
00521 }
00522 //
00523 bool  BoxPartition::setAtomsInOverlapAndCut(){
00524   bool retVal = false; 
00525  // Ajout des atomes des zones overlap et cut :
00526   int nsd, k, nsd_k = 0;
00527   //
00528   double R_overlap= _Rbuffer1, R_cut = _Rbuffer2 ;
00529   //
00530   for(nsd = 0 ; nsd <  _numberOfPartition  ; ++nsd) {
00531     //
00532     //    std::cout << "   Subdomain  " << nsd <<   " Nombre de voisins : "<<_subDomains[nsd].getNumberOfNeighbors()<<std::endl;
00533     for (k = 0 ; k < _subDomains[nsd].getNumberOfNeighbors() ; ++k){ 
00534       nsd_k = _subDomains[nsd].getNumberOfNeighbor(k) ;
00535       if( nsd_k < nsd )
00536         { continue ; }
00537       //      std::cout << "              nsd_k : " << nsd_k << std::endl;
00538       _subDomains[nsd].addAtomsFromSubdomainInOverlapAndCut(k,  _subDomains[nsd_k],*_setOfAtoms, 
00539                                                             R_overlap, R_cut,_useFrag);
00540     }
00541     //for nsd
00542   }
00543   return retVal;
00544 }
00545 //
00546 void BoxPartition::buildLocalNumerotation() {
00547   std::cout << "TO DO : BoxPartition::buildLocalNumerotation(). "<<std::endl;
00548    for(int nsd = 0 ; nsd <  _numberOfPartition  ; ++nsd) {
00549      //     _subDomains[nsd].buildLocalNumerotation(*_setOfAtoms,_atomsInZone) ;
00550    }  //for nsd
00551  }
00552 
00553 
00554 void BoxPartition::sortInFirstDimension(int*  list) {
00555 
00556   // A noter que même si la variable nbAtoms n'a pas besoin d'être un
00557   // argument d'entrée, nous la mettons comme telle, car le partitionneur
00558   // ne restera pas forcément dans la classe CQCAtoms (où la variable
00559   // nbAtoms est connue).
00560 
00561   int      i, j;
00562   int      temp1;
00563   double   temp3;
00564 
00565   for(i = 1 ; i < _setOfAtoms->getNumberOfAtoms() ; ++i) {
00566     temp3 = this->_setOfAtoms->getElement(i)._xyz[0] ;
00567     for(j = 0 ; j < i ; ++j) {
00568 
00569       if (temp3 < (this->_setOfAtoms->getElement(j))._xyz[0] ){
00570         this->_setOfAtoms->swap(i,j) ;
00571         temp1            = list[j];
00572         list[j]          = list[i];
00573         list[i]          = temp1;
00574 #if defined (TT)
00575         temp1            = _charge[j] ;
00576         _charge[j]       = _charge[i] ;  
00577         _charge[i]       = temp1 ;  
00578 #endif
00579       }//if temp3
00580 
00581     }//for j
00582   }//for i
00583 }
00584 bool BoxPartition::computeBoxEquilibrateCharge( const double& Roverlap, const double& RCut) {
00585   //
00586   int     nsd, meanCharge, nb, errorLoad;
00587   double  totalWideSD, totalTmpWideSD, removeLength, gamma, diff;
00588   double* tmpWideSD = new double [_numberOfPartition];
00589   double* wideSD    = new double [_numberOfPartition];
00590   double* c         = new double [_numberOfPartition];
00591   double  RhalfOverlapAndCut = 0.5*Roverlap + RCut ;
00592   double  RhalfOverlap       = 0.5*Roverlap , tmp , error ;
00593   bool isCV ;
00594   //
00595   if( ! (_gridOfPartition[1] == 1 && _gridOfPartition[2] == 1 ) )
00596     { 
00597       std::cerr << " Pas d'equilibrage de charge" << std::endl;
00598       return true; }
00599   isCV        = false ;
00600   //
00601   meanCharge  = 0;
00602   totalWideSD = 0.0;
00603   for(nsd = 0 ; nsd < _numberOfPartition ; ++nsd) {
00604     wideSD[nsd]  = _subDomains[nsd].getXLength() ;
00605     meanCharge  += _subDomains[nsd].load() ;
00606     totalWideSD += wideSD[nsd] ;
00607     c[nsd]       = static_cast<double>(_subDomains[nsd].getNumberOfNeighbors()) ;
00608   }//for i
00609   meanCharge /= _numberOfPartition ; 
00610 
00611   // Calcul des nouveaux wideSD[i] :
00612   // Nouveaux wideSD[i] temporaires :
00613   totalTmpWideSD  = 0; nb = 0 ; removeLength = 0.0 ;
00614   for(nsd = 0 ; nsd < _numberOfPartition ; ++nsd) {
00615     tmp            = c[nsd]*RhalfOverlapAndCut ;
00616     //
00617     // newLi   = newHi - c(i)*( 1/2 Roverlap + Rcut)
00618     //
00619     tmpWideSD[nsd] = (meanCharge / static_cast<double>(_subDomains[nsd].load()))*(wideSD[nsd] + tmp) - tmp;
00620     //    std::cout << "== "<< nsd << " "<< wideSD[nsd] <<" "<< tmpWideSD[nsd] <<std::endl;
00621     //
00622     // Evaluation de la longueur totale temporaire de la molécule :
00623     //
00624     if(tmpWideSD[nsd] <  c[nsd]*RhalfOverlap){
00625       removeLength += c[nsd]*RhalfOverlap ; }
00626     else{
00627       ++nb ; totalTmpWideSD += tmpWideSD[nsd] ;
00628     }
00629   }//for nsd
00630 
00631   // gamma = facteur de réduction :
00632   gamma       = (totalWideSD - removeLength ) / totalTmpWideSD;
00633   totalWideSD = 0.0 ; diff = 0.0; errorLoad = 0 ;
00634   // Les nouveaux wideSD[i] :
00635   for(nsd = 0 ; nsd < _numberOfPartition ; ++nsd) {
00636     if(tmpWideSD[nsd] <  c[nsd]*RhalfOverlap){
00637       tmpWideSD[nsd] =  c[nsd]*RhalfOverlap ;}
00638     else{
00639       tmpWideSD[nsd] = gamma * tmpWideSD[nsd];
00640     }
00641     //    std::cout << "== corection "<< nsd << " "<< wideSD[nsd] <<" "<< tmpWideSD[nsd] <<std::endl;
00642     errorLoad = std::max( errorLoad , meanCharge - _subDomains[nsd].load()) ;
00643     //    std::cout << nsd << " diff load \t" << meanCharge - _subDomains[nsd].load() << "\t  diff dist :  " << fabs(wideSD[nsd] - tmpWideSD[nsd]) << std::endl;
00644     diff  = std::max( diff , fabs(wideSD[nsd] - tmpWideSD[nsd]) );
00645     wideSD[nsd] = tmpWideSD[nsd];
00646     totalWideSD += wideSD[nsd] ;
00647   }//for nsd
00648   error = static_cast<double>(errorLoad) / static_cast<double>(meanCharge) ;
00649   std::cout << std::endl<< std::endl
00650             << "=======================================================================" << std::endl  ;
00651   std::cout << "Equilibrate :   " <<   meanCharge<< std::endl
00652             << "    absolute error on Load     : " << errorLoad << "   error in %   " <<  error*100 << std::endl 
00653             << "    absolute error on distance : " << diff << std::endl << std::endl;
00654   //
00655   //  Convergence si on a une erreur de 5% sur la charge
00656   // 
00657   if( error < 1.0e-2 )
00658     {isCV = true ; 
00659     std::cout << " La methode a convergee." << std::endl;
00660     }
00661   // 
00662   // Move xmin and xmax
00663   //
00664   double s1[3],s2[3] ; 
00665   _subDomains[0].getXmin( s1); 
00666   _subDomains[0].getXmax( s2); 
00667   //
00668   for (int nsd = 0 ; nsd < _numberOfPartition -1; ++nsd){
00669     s1[0] += wideSD[nsd] ;
00670     s2[0] = s1[0] ;
00671     _subDomains[nsd].setXmax( s2); 
00672     _subDomains[nsd+1].setXmin( s1); 
00673     //        std::cout << "=====> New Length : " <<  _subDomains[nsd].getXLength()  << std::endl;
00674   } 
00675   //    std::cout << "=====> New Length : " <<  _subDomains[_numberOfPartition -1].getXLength()  << std::endl;
00676   //  std::cout << "=====>     Total Length  : " <<totalWideSD << std::endl ;
00677   std::cout << "=======================================================================" 
00678               << std::endl<< std::endl;
00679   // Free memory    
00680   delete[] wideSD   ; wideSD = NULL;
00681   delete[] tmpWideSD; tmpWideSD = NULL;
00682   return isCV;
00683 }

Generated on Sat Jan 28 21:07:24 2006 for QC++ by  doxygen 1.4.4