subDomain.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 #include <algorithm>
00019 #include <math.h>
00020 //                                                                                    //
00021 #include "subDomain.h"
00022 #include "commonDef.h"
00023 #include "tools.h"
00024 #include "Atoms.h"
00025 subDomain::subDomain() : _myNum(0),_numberOfAtoms(0), _numberOfAtomsKernel(0),
00026                          _numberOfAtomsOverlap(0),_numberOfAtomsCut(0),
00027                          _indexOfNeighbors(),_indexInKernel(), _atomIsInKernel(),
00028                          _indexInOverlap(NULL),_indexInCut(NULL)
00029 {
00030   _localNum[0]  = _localNum[1]  = _localNum[2]  = 0 ;
00031   _weight[Core_Zone]     = 0;
00032   _weight[Overlap_Zone]  = 0;
00033   _weight[Cut_Zone]      = 0;
00034   //
00035   _xmin[0]          = 0.0 ;
00036   _xmin[1]          = 0.0 ;
00037   _xmin[2]          = 0.0 ;
00038   _xmax[0]          = 0.0 ;
00039   _xmax[1]          = 0.0 ;
00040   _xmax[2]          = 0.0 ;
00041 }
00042 subDomain::~subDomain()
00043 {
00044   if (_indexInOverlap != NULL) {
00045     delete [] _indexInOverlap ; _indexInOverlap = NULL;
00046   }
00047   if (_indexInCut != NULL) {
00048     delete [] _indexInCut ; _indexInCut = NULL;
00049   }
00050 }
00051 void subDomain::resetData(){
00052   _numberOfAtoms = _numberOfAtomsKernel =_numberOfAtomsOverlap = _numberOfAtomsCut = 0 ;
00053   _indexInKernel.clear() ;
00054   _atomIsInKernel.clear() ;
00055   for(unsigned int i = 0 ; i < _indexOfNeighbors.size(); ++i){
00056     _indexInOverlap[i].clear() ;
00057     _indexInCut[i].clear() ;
00058   }
00059   _weight[Core_Zone]     = 0;
00060   _weight[Overlap_Zone]  = 0;
00061   _weight[Cut_Zone]      = 0;
00062   //
00063 }
00064 bool subDomain::initNumberOfOverlapDomain(){
00065   //
00066   bool retVal = true ;
00067   int size    = _indexOfNeighbors.size() ;
00068   if (size == 0){
00069     retVal = false ;}
00070   else {
00071     if (_indexInOverlap == NULL) {
00072       _indexInOverlap = new std::vector<int>[size];
00073     }
00074     else {
00075       retVal = false;
00076     }
00077     if (_indexInCut == NULL) {
00078       _indexInCut     = new std::vector<int>[size];
00079     }
00080     else {
00081       retVal = false;
00082     }
00083   }
00084   return retVal;
00085 }
00086  //
00087 bool subDomain::initNumberOfOverlapDomain(const int& size){
00088   //
00089   bool retVal = true ;
00090   _indexOfNeighbors.resize(size) ;
00091   if (size == 0){
00092     retVal = false ;}
00093   else {
00094     if (_indexInOverlap == NULL) {
00095       _indexInOverlap = new std::vector<int>[size];
00096     }
00097     else {
00098       retVal = false;
00099     }
00100     if (_indexInCut == NULL) {
00101       _indexInCut     = new std::vector<int>[size];
00102     }
00103     else {
00104       retVal = false;
00105     }
00106   }
00107   return retVal;
00108 }
00109  //
00110 //
00111 void subDomain::removeNeigboors(std::vector<int> & list){
00112   //
00113   //  int i, j;
00114 //   for(j= _indexOfNeighbors.size()-1 ; j >=0 ;--j){
00115 //     for(i = 0 ; i< list.size() ){
00116 //       if(_indexOfNeighbors[j] == list[i]){
00117 //      _indexOfNeighbors.erase(v.begin() +j);
00118 //       }
00119 //     }
00120 //   }
00121 }
00122 //
00123 template <class TYPEOFSYSTEM>
00124 bool subDomain::addAtomsFromSubdomainInOverlapAndCut(const int & localOverlapSubDomain, 
00125                                                      subDomain & domain, 
00126                                                      const TYPEOFSYSTEM &atoms, 
00127                                                      //  const setOfAtoms &atoms, 
00128                                                      const double & R_overlap, 
00129                                                      const double & R_cut, bool use_frag) {
00130   int numAtom ;
00131   double bord[3], dist ,  R , R_halfOverlap = 0.5 * R_overlap;
00132   AtomIn  atomNSD1, atomNSD2 ;
00133   int    *localNum;
00134   double coeff[3] ;
00135   //
00136   localNum = domain.getLocalNum() ;
00137   for(int i = 0; i<3 ; ++i){
00138     coeff[i] =  localNum[i] - _localNum[i] ;
00139     if (coeff[i] > 0){
00140       bord[i] = this->_xmax[i];}
00141     else{
00142       bord[i] = this->_xmin[i]; }
00143     //    std::cout <<"  "<<i <<"  " <<  _localNum[i] << "  "<<  localNum[i]<<"  "<< coeff[i]<< "  "<< bord[i]<<std::endl;
00144   }
00145   //  std::cout << "BORD " << bord[0]<< "  "<< bord[1]<< "  "<< bord[2]<< "  "<<std::endl;
00146   R =  0.5* R_overlap + R_cut ;
00147   std::vector<int> & index1           = this->getAtomsIndexInKernel() ;
00148   std::vector<int> & index2           = domain.getAtomsIndexInKernel() ;
00149   std::vector<int> & atomsIsInKernel2 = domain.getAtomsIsInKernel() ;
00150   //
00151   std::vector<int> indexToRemove ;
00152   //
00153   atomNSD1._numDomain = _myNum       ;  
00154   atomNSD1._typeZone  = Overlap_Zone ;
00155   atomNSD1._localNum  = -1 ;
00156   atomNSD2._numDomain = domain. myNumber()  ;  
00157   atomNSD2._typeZone  = Overlap_Zone ;
00158   atomNSD2._localNum  = -1 ;
00159   //
00160   int localOverlapSubDomain1, localOverlapSubDomain2 ;
00161   localOverlapSubDomain1 = localOverlapSubDomain ;
00162   localOverlapSubDomain2 = domain.getLocalNumberOfSundomain( atomNSD1._numDomain);
00163   if(index1.size() == 0 || index2.size() == 0)
00164     {return true;}
00165   //  std::cout <<" After Break "<<index1.size() << "  "<< index2.size()<< std::endl ;
00166   //
00167   // Ajout des atomes de la partie B
00168   //  std::cout << " PARTIE B " << std::endl ;
00169   for(int at = 0 ; at < static_cast<int>(index2.size()) ; ++at){
00170     numAtom = index2[at] ;
00171     //    typename ::atoms::TelementBase & atom  = atoms.getElement(numAtom) ;
00172     //    const atomElement & atom  = atoms.getElement(numAtom) ;
00173     //
00174     dist = fabs((bord[0] - atoms.getElement(numAtom)._xyz[0])*coeff[0]) ;
00175     dist = std::max(dist , fabs((bord[1] - atoms.getElement(numAtom)._xyz[1])*coeff[1])) ;
00176     dist = std::max(dist , fabs((bord[2] - atoms.getElement(numAtom)._xyz[2])*coeff[2])) ;
00177     //
00178     if ( dist > R )
00179       { continue ;}
00180     //    std::cout <<at <<" Atome "<< numAtom << "  " << dist << "   " << R <<std::endl;
00181     if ( dist < R_halfOverlap ){
00182       atomsIsInKernel2[at] = 0 ;
00183       indexToRemove.push_back(at) ;
00184       this->addIndexInOverlap(localOverlapSubDomain1,numAtom) ;
00185       domain.addIndexInOverlap(localOverlapSubDomain2,numAtom) ;
00186     }
00187     else {
00188       atomNSD1._typeZone = Cut_Zone ;
00189       this->addIndexInCut(localOverlapSubDomain1,numAtom)  ;
00190       atomNSD1._typeZone = Overlap_Zone ;
00191     }
00192   }
00193   std::sort(indexToRemove.begin(),indexToRemove.end());
00194 //    std::cout <<"INDEX 2 : "  ;
00195 //    for(unsigned int ii = 0 ; ii < index2.size() ; ++ii ){
00196 //      std::cout << index2[ii] << "  " ;
00197 //    }
00198 //    std::cout << std::endl ;
00199 //    std::cout <<"INDEX 2 TO REMOVE :  "  ;
00200 //    for(unsigned int ii = 0 ; ii < indexToRemove.size() ; ++ii ){
00201 //      std::cout << indexToRemove[ii] << "  " ;
00202 //    }
00203 //    std::cout << std::endl ;
00204  // CAS LINEAIRE
00205   std::vector<int>::iterator startIterator = index2.begin(); 
00206   index2.erase(startIterator,startIterator + indexToRemove.size());
00207  //
00208 //   std::cout <<"INDEX 2 : "  ;
00209 //   for(unsigned int ii = 0 ; ii < index2.size() ; ++ii ){
00210 //     std::cout << index2[ii] << "  " ;
00211 //   }
00212 //   std::cout << std::endl ;
00213   indexToRemove.clear() ;
00214   //
00215   //
00216   // Ajout des atomes de la partie A 
00217   // std::cout << std::endl<<std::endl<<" PARTIE A " << std::endl ;
00218   for(int at = 0 ; at < static_cast<int>(index1.size()) ; ++at){
00219    numAtom = index1[at] ;
00220    //   const atomElement & atom  = atoms.getElement(numAtom) ;
00221    //
00222     dist = fabs( (bord[0] - atoms.getElement(numAtom)._xyz[0])*coeff[0] ) ;
00223     dist = std::max( dist , fabs( (bord[1] - atoms.getElement(numAtom)._xyz[1])*coeff[1]) ) ;
00224     dist = std::max( dist , fabs( (bord[2] - atoms.getElement(numAtom)._xyz[2])*coeff[2] )) ;
00225 //     dist = fabs( (bord[0] - atom._xyz[0])*coeff[0] ) ;
00226 //     dist = std::max( dist , fabs( (bord[1] - atom._xyz[1])*coeff[1]) ) ;
00227 //     dist = std::max( dist , fabs( (bord[2] - atom._xyz[2])*coeff[2] )) ;
00228     // 
00229     if ( dist > R )
00230       { continue ;}
00231     //    std::cout <<at <<" Atome "<< numAtom << "  " << dist << "   " << R <<std::endl;
00232     if ( dist < R_halfOverlap ){
00233       this->addIndexInOverlap(localOverlapSubDomain1,numAtom) ;
00234       _atomIsInKernel[at] = 0; // on le marque pour l'enlever du kernel
00235       indexToRemove.push_back(at) ;
00236       domain.addIndexInOverlap(localOverlapSubDomain2,numAtom) ;
00237     }
00238     else {
00239       atomNSD2._typeZone = Cut_Zone ;
00240       domain.addIndexInCut(localOverlapSubDomain2,numAtom)  ;
00241       atomNSD2._typeZone = Overlap_Zone ;
00242     }
00243   }
00244   //
00245   // Enlever les atomes du kernel qui sont mis dans l'overlap
00246   //
00247 //   std::sort(indexToRemove.begin(),indexToRemove.end());
00248 //   std::cout <<"INDEX 1 : "  ;
00249 //   for(unsigned int ii = 0 ; ii < index1.size() ; ++ii ){
00250 //     std::cout << index1[ii] << "  " ;
00251 //   }
00252 //   std::cout << std::endl ;
00253 //   std::cout <<"INDEX 1 TO REMOVE : "  ;
00254 //   for(unsigned int ii = 0 ; ii < indexToRemove.size() ; ++ii ){
00255 //     std::cout << indexToRemove[ii] << "  " ;
00256 //   }
00257 //   std::cout << std::endl ;
00258  // CAS LINEAIRE
00259  startIterator = index1.end()  - indexToRemove.size(); 
00260  index1.erase(startIterator, index1.end());
00261  //
00262 //   std::cout <<"INDEX 1 : "  ;
00263 //   for(unsigned int ii = 0 ; ii < index1.size() ; ++ii ){
00264 //     std::cout << index1[ii] << "  " ;
00265 //   }
00266 //   std::cout << std::endl ;
00267  indexToRemove.clear() ;  
00268  return true ;
00269 }
00270 void subDomain::computeNumberOfAtoms(){
00271   //
00272   unsigned int i ;
00273   _numberOfAtomsKernel =_numberOfAtomsOverlap =  _numberOfAtomsCut = 0 ;
00274   for( i = 0 ; i < _indexOfNeighbors.size() ; ++i){
00275     _numberOfAtomsOverlap += _indexInOverlap[i].size() ;
00276     _numberOfAtomsCut     += _indexInCut[i].size() ;
00277   }
00278   _numberOfAtomsKernel = _indexInKernel.size() ;
00279 
00280   _numberOfAtoms = _numberOfAtomsOverlap + _numberOfAtomsCut + _numberOfAtomsKernel;
00281 
00282   std::cout << _myNum << " numberOfAtomsKernel : "<< _numberOfAtomsKernel
00283             <<"  numberOfAtomsOverlap = "<< _numberOfAtomsOverlap
00284             << "  numberOfAtomsCut =  "<< _numberOfAtomsCut
00285             <<" numberOfAtoms = "<<  _numberOfAtoms<< std::endl;
00286 }
00287 
00288 inline void  setLocalNum(std::vector<AtomIn> & vec ,const AtomIn &atom){
00289   for(unsigned int i = 0 ; i < vec.size() ; ++i){
00290     if(vec[i]._numDomain == atom._numDomain && (vec[i]._typeZone == atom._typeZone)){
00291       vec[i]._localNum = atom._localNum;    
00292     }
00293   }
00294 }
00295 template <class TYPEOFSYSTEM>
00296 void subDomain::buildLocalNumerotation(const  TYPEOFSYSTEM &atoms, std::vector<AtomIn> *atomsInZone) 
00297 {
00298   int localNum, atK, atO, atOIn[_indexOfNeighbors.size()];
00299   int  i, imin=0, j ;
00300   double x=0., xInOver[_indexOfNeighbors.size()];
00301   bool testInOver ;
00302   //
00303   localNum = 0 ;
00304   //
00305   //  std::vector<int> atomInKernel( _numberOfAtomsKernel) ;
00306   //
00307   atK = atO               = 0 ;
00308   for( i = 0 ; i < static_cast<int>(_indexOfNeighbors.size()) ; ++i ){
00309     std::sort(_indexInOverlap[i].begin(), _indexInOverlap[i].end());
00310     atOIn[i]   = 0 ;
00311     xInOver[i] = atoms.getElement(_indexInOverlap[i][atOIn[i]])._xyz[0] ;
00312     //    std::cout << i << "  " <<  atOIn[i] << "  "<<  xInOver[i] << std::endl;
00313   }
00314   //  std::cout << std::endl;
00315   AtomIn atomK, atomO ;
00316   //
00317   atomK._numDomain =  atomO._numDomain = _myNum;
00318   atomK._typeZone  =  Core_Zone ;   atomO._typeZone  =  Overlap_Zone ; 
00319   //
00320   testInOver = true ; 
00321   bool stopK,stopO ;
00322   //
00323   stopK = stopO = false ;
00324   // 
00325   //  std::cout <<" _atomIsInKernel[atK ] == 0 " << _atomIsInKernel[atK ] <<std::endl;
00326   //
00327   //  while(atK <static_cast<int>(_indexInKernel.size()) && _atomIsInKernel[atK ] == 0   )
00328   while(atK <static_cast<int>(_indexInKernel.size()) )
00329     {atK++; }
00330   if(atK < static_cast<int>(_indexInKernel.size()) ){  
00331     x    =  atoms.getElement(_indexInKernel[atK])._xyz[0] ;
00332   }
00333   //
00334   //  while ( !(atK  <  _indexInKernel.size()) && !( atO <  _numberOfAtomsOverlap) ) {
00335   while ( ( localNum <  _numberOfAtoms) ) {
00336     //
00337     // find the minimum in  overlap zones
00338     if(testInOver){
00339       imin = 0 ; 
00340       for( i = 1 ; i < static_cast<int>(_indexOfNeighbors.size()) ; ++i ){
00341         if( xInOver[imin] > xInOver[i])
00342           { imin = i ; }
00343       }
00344       testInOver = false ;
00345     } 
00346     //    std::cout << localNum << "   atK " << atK << " x  "<< x << "  atO " <<atO <<"  " << imin<< "  "<< xInOver[imin]<<std::endl;
00347     if (xInOver[imin] > x ){
00348       // l'atome dans le kernel est plus petit 
00349       //       std::cout << "l'atome dans le kernel est plus petit " << std::endl ;
00350       atomK._localNum  =  localNum ; 
00351       setLocalNum(atomsInZone[_indexInKernel[atK]],atomK) ; 
00352       localNum++;  atK++;
00353       x     = 999999999.0;
00354       //      while(atK <static_cast<int>(_indexInKernel.size()) && _atomIsInKernel[atK ] == 0   )
00355       while(atK <static_cast<int>(_indexInKernel.size()) )
00356         {atK++; }
00357       
00358       if(atK <static_cast<int>(_indexInKernel.size()) ){  
00359         x    =  atoms.getElement(_indexInKernel[atK])._xyz[0] ;
00360       }
00361       else{
00362         stopK = true;}
00363     }
00364     else { 
00365       // l'atome dans l'Overlap est plus petit
00366       //
00367       atomO._localNum  =  localNum ; 
00368       setLocalNum(atomsInZone[(_indexInOverlap[imin])[atOIn[imin]]],atomO ); 
00369       localNum++; 
00370       // deplacement des curseurs
00371       atO++; atOIn[imin]++;   
00372       if(atOIn[imin]< static_cast<int>(_indexInOverlap[imin].size()) ) {
00373         xInOver[imin] = atoms.getElement(_indexInOverlap[imin][atOIn[imin]])._xyz[0] ;}
00374       else{ 
00375         xInOver[imin] = 99999999.0 ; 
00376         if( atO == _numberOfAtomsOverlap )
00377           {stopO = true;}}
00378       testInOver = true ;
00379     }
00380     if( stopK && stopO)
00381       {break ;}
00382   }
00383   //
00384   //  Numerotation des zone de cut
00385   //
00386   atomO._typeZone  =  Cut_Zone ; atomO._localNum  =  localNum - 1 ;;
00387   for( i = 0 ; i < static_cast<int>(_indexOfNeighbors.size()) ; ++i ){
00388     for( j = 0 ; j < static_cast<int>(_indexInCut[i].size()) ; ++j){
00389       atomO._localNum++ ;
00390       //      std::cout << atomO._localNum << "   atC " << _indexInCut[i][j]<<std::endl;
00391       setLocalNum(atomsInZone[(_indexInCut[i])[j]],atomO) ; 
00392     } 
00393   }
00394   //    std::cout << " local Num " <<  atomO._localNum + 1  << " _numberOfAtoms " <<_numberOfAtoms <<std::endl;
00395 }
00396 
00397 // void subDomain::writeINXYZ(ofstream file, const setOfAtoms &atoms){
00398 
00399 
00400 // }
00401 
00402 template bool subDomain::addAtomsFromSubdomainInOverlapAndCut< setOfAtoms >(const int & localOverlapSubDomain, 
00403                                                                   subDomain & domain, 
00404                                                                   const setOfAtoms &atoms, 
00405                                                                   const double & R_overlap, const double & R_cut, bool use_frag) ;
00406 
00407 template void subDomain::buildLocalNumerotation< setOfAtoms >(const  setOfAtoms &atoms, std::vector<AtomIn> *atomsInZone) ; 

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