AApartition.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 "AApartition.h"
00020 #include "tools.h"
00021 //
00022 #include <iomanip>
00023 #include <vector>
00024 #include <algorithm>
00025 //
00026 AAPartition::AAPartition() : PartitionBase< PDBstructure >()
00027 {
00028 }
00029 
00030 AAPartition::AAPartition(const int size, const double R1, const double R2) : PartitionBase< PDBstructure >(size,R1,R2)
00031 {
00032 }
00033 
00034 AAPartition::~AAPartition()
00035 {}
00036 
00037 bool  AAPartition::partitionate(const bool use_frag){
00038   //
00039   _useFrag =  use_frag ;
00040   //
00041   std::cout << "partitionate "<< std::endl;
00042   double coinMin[3] , coinMax[3], LengthOfMolecule[3] ;
00043   std::cout << "   _setOfAtoms->findMinMax "<< std::endl;
00044   this->_setOfAtoms->findMinMax(coinMin , coinMax  ) ;
00045   //
00046   LengthOfMolecule[0]  =   coinMax[0] - coinMin[0] ;
00047   LengthOfMolecule[1]  =   coinMax[1] - coinMin[1] ;
00048   LengthOfMolecule[2]  =   coinMax[2] - coinMin[2] ;
00049   //
00050   //
00051   std::cout << " XMIN            : "<< coinMin[0] << "  "<< coinMin[1] << "  "<< coinMin[2] << std::endl;
00052   std::cout << " XMAX            : "<< coinMax[0] << "  "<< coinMax[1] << "  "<< coinMax[2] << std::endl;
00053   std::cout << "LengthOfMolecule : " <<LengthOfMolecule[0]<< "  "<<LengthOfMolecule[1]
00054             << "  "<<LengthOfMolecule[2]<< "  "<<std::endl;
00055   //
00056   std::cout << "  ---> setNeigboors " << std::endl;
00057   //
00058   for( int nsd = 0 ; nsd < _numberOfPartition  ; ++nsd){
00059     _subDomains[nsd].setMyNumber(nsd);
00060   }
00061   //
00062   // Add atoms in the kernel
00063   //
00064   std::cout << "  ---> set Atoms In Kernel " << std::endl;
00065   this->setAtomsInKernel() ;
00066   //
00067   this->setNeigboors();
00068   //
00069   //
00070   if(use_frag){
00071     build_fragmentStructure( ) ;
00072   }
00073   //
00074   // Add atoms in Overlap and Cut
00075   //
00076   std::cout << "  ---> set Atoms In Overlap and Cut " << _Rbuffer1 << "  "<<  _Rbuffer2 <<std::endl;
00077   this->setAtomsInOverlapAndCut(); 
00078   //
00079   std::cout << "  ---> set number Of Atoms in subdomains  " <<std::endl;
00080    for (int nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00081      _subDomains[nsd].computeNumberOfAtoms();
00082    }
00083   std::cout << "  ---> set Load of subdomains "  <<std::endl;
00084   this->setLoad() ;
00085   //
00086   std::cout << "  ---> buildLocalNumerotation "  <<std::endl;
00087   //
00088   this->buildLocalNumerotation(); // BUG ????
00089   //
00090   std::cout << "  ---> statistiques "  <<std::endl;
00091   this->statistiques() ;
00092   //
00093   // For QCPP
00094   //
00095   this->buildAtomsInZone() ;
00096   //
00097   //
00098   return true ;
00099 }
00100 
00101 void AAPartition::setNeigboors(){
00102   //
00103   double  atomMax[3], atomMin[3];
00104   for( int nsd = 0 ; nsd <  _numberOfPartition  ; ++nsd) {
00105     //
00106     atomMax[0] = atomMax[1] = atomMax[2] = -10000.0 ;
00107     atomMin[0] = atomMin[1] = atomMin[2] =  10000.0 ;
00108     std::vector<int>& index      = _subDomains[nsd].getAtomsIndexInKernel();
00109     //
00110     int i ;
00111     for ( unsigned int j = 0 ; j < index.size() ;  ++j){
00112       i          =  index[j] ;
00113       atomMax[0] = std::max(atomMax[0], _setOfAtoms->getElement(i)._xyz[0]) ;
00114       atomMax[1] = std::max(atomMax[1], _setOfAtoms->getElement(i)._xyz[1]) ;
00115       atomMax[2] = std::max(atomMax[2], _setOfAtoms->getElement(i)._xyz[2]) ;
00116       atomMin[0] = std::min(atomMin[0], _setOfAtoms->getElement(i)._xyz[0]) ;
00117       atomMin[1] = std::min(atomMin[1], _setOfAtoms->getElement(i)._xyz[1]) ;
00118       atomMin[2] = std::min(atomMin[2], _setOfAtoms->getElement(i)._xyz[2]) ;
00119     }
00120     _subDomains[nsd].setXmin(atomMin);
00121     _subDomains[nsd].setXmax(atomMax);
00122   } //for nsd
00123   //
00124   //
00125   double dist1,dist2,dist ; 
00126 //   int ipos[3] ;
00127 //   int jpos[3] ;
00128   double C1[3], C2[3], R =  _Rbuffer1 + _Rbuffer2;
00129   //
00130   C1[0] = C1[1] = C1[2] = C2[0] = C2[1] = C2[2] = 0.0;
00131   for(int nsd = 0 ; nsd <  _numberOfPartition  ; ++nsd) {
00132     //  atomMin = centre nsd
00133     //  atomMax = centre nsd1
00134     _subDomains[nsd].getCentre(C1) ;
00135     _subDomains[nsd]. getXmin(atomMin) ;
00136     _subDomains[nsd]. getXmax(atomMax) ;
00137 //     std::cout << "Sous domaine : "<< nsd << std::endl;
00138 //     std::cout << " XMIN            : "<< atomMin[0] << "  "<< atomMin[1] << "  "<< atomMin[2] << std::endl;
00139 //     std::cout << " XMAX            : "<< atomMax[0] << "  "<< atomMax[1] << "  "<< atomMax[2] << std::endl;
00140 //     std::cout << " CENTRE          : "<<C1[0]<< "  "<<C1[1]      << "  "<<C1[2]<< "  "<<std::endl;
00141      dist1 =  distance(atomMin,C1) ;
00142 //     std::cout << "                         dist1 " <<dist1<<std::endl;
00143     for(int nsd1 = 0 ; nsd1 <  nsd  ; ++nsd1) {
00144       _subDomains[nsd1].getCentre(C2) ;
00145       _subDomains[nsd1]. getXmin(atomMin) ;
00146       dist  = distance(C1,C2) ;
00147       dist2 = distance(atomMin,C2) ;
00148       //    std::cout << dist1 << "  "<< dist << "  " << dist2 << " "<< dist - dist1 - dist2 -R<<std::endl;
00149       if(dist - dist1 - dist2 < R ){
00150         //              std::cout << "nsd : "<< nsd << " et nsd& "<< nsd1 <<" sont voisins " << std::endl;
00151         _subDomains[nsd].addNeigboors(nsd1);
00152         _subDomains[nsd1].addNeigboors(nsd);
00153       }
00154     }// nsd1
00155   } //nsd 
00156   
00157   //
00158   //  Allocation
00159   //
00160   for( int nsd = 0 ; nsd < _numberOfPartition  ; ++nsd){
00161     _subDomains[nsd].initNumberOfOverlapDomain() ;
00162     //    std::cout << "SD "<< nsd<< " a " <<  _subDomains[nsd].getNumberOfNeighbors() << "  voisins. "
00163     //        << _subDomains[nsd].getNumberOfNeighbors() << std::endl; 
00164   }
00165 }
00166 //
00167  bool  AAPartition::setAtomsInKernel(){
00168 
00169  bool retVal = false;
00170  int numberOfAtoms = this->_setOfAtoms->numberOfAtoms(), nsd, nbAtoms ;
00171  //
00172  nsd = nbAtoms = 0 ;
00173  while(nbAtoms < numberOfAtoms){
00174    nsd =  this->_setOfAtoms->getElement(nbAtoms).numberOfResidu() ;
00175    if(nsd < _numberOfPartition){
00176      _subDomains[nsd].getAtomsIndexInKernel().push_back( nbAtoms );
00177      _AAtoSD[nbAtoms] = nsd ;
00178    }
00179    else{
00180      std::cerr << "Error current subdomain "<<nsd <<" > number of subdomains "<<  _numberOfPartition <<std::endl;
00181      exit(EXIT_FAILURE) ;
00182  }
00183    nbAtoms++;
00184  }
00185  //
00186  return retVal ;
00187 }
00188 //
00189 // bool  AAPartition::setAtomsInOverlapAndCut(){
00190 //   bool retVal = false; 
00191 //   // Ajout des atomes des zones overlap et cut :
00192 //   int k, nsd, nsd_k = 0, atomI, atomJ ;
00193 //   double dist=0., distMin ;
00194 //   double R = _Rbuffer1 + _Rbuffer2, R1 =  _Rbuffer1* _Rbuffer1 ;
00195 //   //
00196 //   R = R*R ;
00197 //   for(nsd = 0 ; nsd <  _numberOfPartition  ; ++nsd) {
00198 //     //
00199 //     std::vector<int>& index      = _subDomains[nsd].getAtomsIndexInKernel();
00200 //     for (k = 0 ; k < _subDomains[nsd].getNumberOfNeighbors() ; ++k){ 
00201 //       nsd_k = _subDomains[nsd].getNumberOfNeighbor(k) ;
00202 //       //      std::cout << "              nsd_k : " << nsd_k << std::endl;
00203 //       std::vector<int>& indexZONE1 = _subDomains[nsd].getAtomsIndexInOverlap(k);
00204 //       std::vector<int>& indexZONE2 = _subDomains[nsd].getAtomsIndexInCut(k);
00205 //       //
00206 //       std::vector<int>& index_k      = _subDomains[nsd_k].getAtomsIndexInKernel();
00207 //       //
00208 //       // Boucle sur les atomes
00209 //       //
00210 //       for ( unsigned int j = 0 ; j < index_k.size() ; ++j){
00211 //      atomJ               = index_k[j] ;
00212 //      //      atomElement & atome = this->_setOfAtoms->getElement(atomJ) ;
00213 //      distMin             = 9999. ;
00214 //      for ( unsigned int i = 0 ; i < index.size() ; ++i){
00215 //        dist    = this->_setOfAtoms->getElement(atomJ).dist2To(this->_setOfAtoms->getElement(index[i]).pos() );
00216 //        if (distMin > dist){
00217 //          distMin = dist ; atomI = index[i] ;
00218 //        }
00219 //      } // for i
00220 //      if( distMin > R){
00221 //        continue ; }
00222 //      else {
00223 //        if (distMin < R1 ){
00224 //          //      indexZONE1.push_back( this->_setOfAtoms->getElement(atomJ).number()) ;
00225 //          indexZONE1.push_back(atomJ) ;
00226 //        }
00227 //        else {
00228 //          //      indexZONE2.push_back( this->_setOfAtoms->getElement(atomJ).number()) ;
00229 //          indexZONE2.push_back(atomJ) ;
00230 //        }
00231 //      }
00232 //       } // for j
00233 //     }  //for k
00234 //   } //for nsd
00235 //   return retVal;
00236 // }
00237 
00238 //
00239 void AAPartition::buildLocalNumerotation() {
00240   for(int nsd = 0 ; nsd <  _numberOfPartition  ; ++nsd) {
00241     //    _subDomains[nsd].buildLocalNumerotation(*_setOfAtoms,_atomsInZone) ;
00242   }  //for nsd
00243 }
00244 //
00245 void AAPartition::sortInFirstDimension(int*  list) {
00246 
00247   // A noter que même si la variable nbAtoms n'a pas besoin d'être un
00248   // argument d'entrée, nous la mettons comme telle, car le partitionneur
00249   // ne restera pas toujours dans la classe CQCAtoms (où la variable
00250   // nbAtoms est connue).
00251 
00252   int      i, j;
00253   int      temp1;
00254   double   temp3;
00255 
00256   for(i = 1 ; i < _setOfAtoms->numberOfAtoms() ; ++i) {
00257     temp3 = this->_setOfAtoms->getElement(i)._xyz[0] ;
00258     for(j = 0 ; j < i ; ++j) {
00259 
00260       if (temp3 < (this->_setOfAtoms->getElement(j))._xyz[0] ){
00261         this->_setOfAtoms->swap(i,j) ;
00262         temp1            = list[j];
00263         list[j]          = list[i];
00264         list[i]          = temp1;
00265       }//if temp3
00266 
00267     }//for j
00268   }//for i
00269 }
00270 bool AAPartition::computeAAEquilibrateCharge( const double& Roverlap, const double& RCut) {
00271   //
00272   return true;
00273 }
00274 
00275 bool  AAPartition::valid(){
00276   return false;
00277 }
00278 //
00279 void  AAPartition::statistiques(){
00280   //
00281   // Charge moyenne
00282   //
00283   int meanLoad = 0 ;
00284   double ecart = 0.0 ;
00285   for( int nsd = 0 ; nsd < _numberOfPartition  ; ++nsd){
00286     meanLoad +=  _subDomains[nsd].load() ;
00287   }
00288   meanLoad /= _numberOfPartition ;
00289   for( int nsd = 0 ; nsd < _numberOfPartition  ; ++nsd){
00290     ecart +=  (_subDomains[nsd].load() - meanLoad)*(_subDomains[nsd].load() - meanLoad) ;
00291   }
00292   ecart = std::sqrt(ecart) ;
00293   //
00294   // on affiche sous la forme d'un tableau
00295   //
00296   std::cout <<std::endl<< std::endl;
00297   std::cout <<"Number Of Subdomains : "<<_numberOfPartition<< "  Average Load : " << meanLoad
00298             <<" ecart type : " << ecart<< std::endl ;
00299   // 
00300   //  id residu chargeTotale Charge{ Kernel , B1 , B2 } ecartType
00301   for( int nsd = 0 ; nsd < _numberOfPartition  ; ++nsd){
00302     
00303     std::cout << nsd << "  " 
00304               << this->_setOfAtoms->getElement(_subDomains[nsd].getAtomsIndexInKernel()[0]).getResiduName() 
00305               << "  "  <<  _subDomains[nsd].load() <<"  "<<  _subDomains[nsd].load() -meanLoad<<"  "
00306               << _subDomains[nsd].loadInKernel()<<"  "
00307               << _subDomains[nsd].loadInOverlap()<<"  "<< _subDomains[nsd].loadInCut() << std::endl ;
00308   }
00309   //
00310 
00311   //
00312 }
00313 //

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