00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <stdlib.h>
00019 #include <iomanip>
00020 #include <iostream>
00021
00022
00023 #include "PartitionBaseT.h"
00024
00025 #include "tools.h"
00026 #include "Atoms.h"
00027 #include "pdbstruct.h"
00028
00029
00030 template <class TYPEOFATOMS>
00031 PartitionBase<TYPEOFATOMS>::PartitionBase() : _partitionType("QC_Standard_Part"),_numberOfPartition(0),_charge(NULL),_setOfAtoms(NULL),
00032 _setOfFragments(NULL),
00033 _subDomains(NULL),_atomsInZone(NULL),_Rbuffer1(3.0),_Rbuffer2(1.5),
00034 _useFrag(false)
00035 {}
00036 template <class TYPEOFATOMS>
00037 PartitionBase<TYPEOFATOMS>::PartitionBase(const int size, const double R1, const double R2) :
00038 _partitionType("QC_Standard_Part"),_numberOfPartition(size),_charge(NULL),_setOfAtoms(NULL),
00039 _setOfFragments(NULL),_subDomains(new subDomain[size]),
00040 _atomsInZone(NULL),_Rbuffer1(R1),_Rbuffer2(R2),
00041 _useFrag(false)
00042 {}
00043
00044 template <class TYPEOFATOMS>
00045 PartitionBase<TYPEOFATOMS>::~PartitionBase()
00046 {
00047 _setOfAtoms = NULL;
00048 if(_subDomains != NULL)
00049 {delete [] _subDomains ; _numberOfPartition = 0 ;}
00050 if(_charge != NULL)
00051 {delete [] _charge ; _charge = NULL ;}
00052 }
00053
00054 template <class TYPEOFATOMS>
00055 void PartitionBase<TYPEOFATOMS>::resetData(){
00056
00057 std::cout << " ---> resetData " <<std::endl;
00058 for (int nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00059 _subDomains[nsd].resetData();
00060 }
00061
00062 delete [] _atomsInZone ;
00063 _atomsInZone = new std::vector<AtomIn>[this->_setOfAtoms->numberOfAtoms()];
00064 }
00065 template <class TYPEOFATOMS>
00066 bool PartitionBase<TYPEOFATOMS>::attachAtoms( TYPEOFATOMS *atoms )
00067 {
00068 if( atoms == NULL){
00069 std::cerr <<"PartitionBase::attachAtoms :No Atoms to attach. " <<std::endl ;
00070 exit(1) ;
00071 }
00072 _useFrag = false ;
00073 _setOfAtoms = atoms ;
00074 int nbAtoms = _setOfAtoms ->numberOfAtoms() ;
00075
00076 _charge = new int[_setOfAtoms ->numberOfAtoms() ] ;
00077 int num ;
00078 for ( int i = 0 ; i< _setOfAtoms ->numberOfAtoms() ; ++i){
00079 num = giveNumberOfOA(_setOfAtoms->getAtomName(i)) ;
00080 if ( num != 0) {
00081 _charge[i] = num ;
00082 }
00083
00084
00085
00086
00087
00088 }
00089 _atomsInZone = new std::vector<AtomIn>[nbAtoms];
00090
00091 return true ;
00092 }
00093 template <class TYPEOFATOMS>
00094 bool PartitionBase<TYPEOFATOMS>::attachAtoms( setOfFragments *fragments )
00095 {
00096 if( fragments == NULL){
00097 std::cerr <<"PartitionBase::attachAtoms :No Atoms to attach. " <<std::endl ;
00098 exit(1) ;
00099 }
00100 std::cerr <<"PartitionBase::attachAtoms :No Atoms to attach. " <<std::endl ;
00101 exit(1) ;
00102 return true ;
00103 }
00104 template <class TYPEOFATOMS>
00105 bool PartitionBase<TYPEOFATOMS>::build_fragmentStructure( ){
00106
00107 _nextATofFrag.resize(this->_setOfAtoms->numberOfAtoms()) ;
00108 _firstATofFrag.resize(this->_setOfAtoms->numberOfFragments()) ;
00109
00110 for(unsigned int i = 0 ; i < _nextATofFrag.size() ; ++i){
00111 _nextATofFrag[i] = -1 ;
00112 }
00113 for(unsigned int i = 0 ; i < _firstATofFrag.size() ; ++i){
00114 _firstATofFrag[i] = -1 ;
00115 }
00116
00117
00118
00119 int fragNum, j ;
00120 for (int i = 0 ; i < this->_setOfAtoms->numberOfAtoms() ; ++i){
00121 fragNum = _setOfAtoms->getElement(i).getFragmentNumber() -1;
00122
00123 if( _firstATofFrag[fragNum] == -1){
00124 _firstATofFrag[fragNum] = i ;}
00125 else{
00126 j = _firstATofFrag[fragNum] ;
00127
00128 while ( _nextATofFrag[j] != -1){
00129 j = _nextATofFrag[j] ;
00130
00131 }
00132
00133 _nextATofFrag[j] = i ;
00134 }
00135 }
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 return true;
00149 }
00150
00151
00152 template <class TYPEOFATOMS>
00153 bool PartitionBase<TYPEOFATOMS>::fragmentsToAtoms()
00154 {
00155 std::cout << "BUGS BUGS _atomsInZone[] not allocated " <<std::endl;
00156 if(!_useFrag)
00157 {return true ;}
00158
00159 int nbAtoms = _setOfFragments->totalNumberOfAtoms() ;
00160
00161 std::vector<AtomIn> * atomsInZone ;
00162 atomsInZone = new std::vector<AtomIn>[nbAtoms];
00163
00164 int *posInAtoms, *currentLocalNumber;
00165 _setOfAtoms = new TYPEOFATOMS(nbAtoms) ;
00166 posInAtoms = new int[_setOfFragments-> numberOfFraments() ] ;
00167 currentLocalNumber = new int[_numberOfPartition ];
00168
00169 int i, j, nsd, numberOfAtomsInFragment, currentAtoms;
00170 unsigned int k ;
00171 posInAtoms[0] = 0 ;
00172 for ( i = 1 ; i< _setOfFragments->numberOfFraments() ; ++i){
00173 posInAtoms[i] = posInAtoms[0] + _setOfFragments->getElement(i-1).numberOfAtoms();
00174 }
00175
00176 for ( nsd = 1 ; nsd< _numberOfPartition ; ++nsd){
00177 currentLocalNumber[nsd] = 0;
00178 }
00179
00180
00181
00182 currentAtoms = 0 ;
00183
00184 for ( i = 0 ; i< _setOfFragments->numberOfFraments() ; ++i){
00185 FragmentElement & frag = _setOfFragments->getElement(i);
00186 numberOfAtomsInFragment = frag.numberOfAtoms() ;
00187
00188
00189 for (j = 0; j < numberOfAtomsInFragment ; ++j){
00190
00191 for( k = 0 ; k < _atomsInZone[j].size() ; ++k){
00192 atomsInZone[currentAtoms].push_back(_atomsInZone[i][k]) ;
00193 }
00194 ++currentAtoms;
00195 }
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 }
00212
00213 return true ;
00214 }
00215
00216 template <class TYPEOFATOMS>
00217 void PartitionBase<TYPEOFATOMS>::setLoad(){
00218
00219 int nsd, k, loadOverlap, loadCut,loadKernel;
00220 unsigned int at;
00221
00222
00223 for(nsd = 0 ; nsd < _numberOfPartition ; ++nsd) {
00224
00225 loadOverlap = loadCut = loadKernel = 0 ;
00226 std::vector<int> & indexK = _subDomains[nsd].getAtomsIndexInKernel() ;
00227
00228 for (k = 0 ; k < _subDomains[nsd].getNumberOfNeighbors() ; ++k){
00229 std::vector<int> & indexO = _subDomains[nsd].getAtomsIndexInOverlap(k) ;
00230 std::vector<int> & indexC = _subDomains[nsd].getAtomsIndexInCut(k) ;
00231 for(at = 0 ; at < indexO.size() ; ++at){
00232 loadOverlap += _charge[indexO[at]] ;
00233 }
00234 for(at = 0 ; at < indexC.size() ; ++at){
00235 loadCut += _charge[indexC[at]] ;
00236 }
00237 }
00238 for(at = 0 ; at < indexK.size() ; ++at){
00239 loadKernel += _charge[indexK[at]] ;
00240 }
00241
00242
00243
00244
00245
00246
00247 _subDomains[nsd].loadInKernel() = loadKernel ;
00248 _subDomains[nsd].loadInOverlap() = loadOverlap ;
00249 _subDomains[nsd].loadInCut() = loadCut ;
00250 }
00251 }
00252
00253 template <class TYPEOFATOMS>
00254 void PartitionBase<TYPEOFATOMS>::fixAtomsInZone(){
00255
00256 bool toErase ;
00257 for( int i = 0 ; i< this->_setOfAtoms->numberOfAtoms() ; ++i ){
00258 if(_atomsInZone[i].size() <= 2){
00259 continue ;}
00260 toErase = false ;
00261 if((_atomsInZone[i])[0]._localNum == -1 ){
00262
00263 toErase = true ;
00264 }
00265 if(toErase)
00266 {(_atomsInZone[i]).erase(_atomsInZone[i].begin() ) ;}
00267 }
00268 }
00269
00270 template <class TYPEOFATOMS>
00271 void PartitionBase<TYPEOFATOMS>::writeInFile(std::ofstream &out)
00272 {
00273
00274
00275 int i , nsd ;
00276
00277 int j ;
00278 std::vector< std::vector<int> > numOfSubDom( _numberOfPartition);
00279 std::vector< int > load(_numberOfPartition);
00280 int loadAtoms = 0 , numSubDom;
00281 for(i=0; i < _numberOfPartition ; ++i){
00282 numOfSubDom[i].resize( _numberOfPartition) ;
00283 load[i] = 0 ;
00284 for(j = 0; j < _numberOfPartition ; ++j){
00285 (numOfSubDom[i])[j] = 0 ;
00286 }
00287 }
00288 int numSubDom1 ;
00289 for(i = 0 ; i < _setOfAtoms->numberOfAtoms() ; ++i) {
00290
00291 loadAtoms = giveNumberOfOA(_setOfAtoms->getAtomName(i)) ;
00292 std::vector<AtomIn> & v = _atomsInZone[i] ;
00293 for (unsigned int k = 0 ; k < v.size() ; ++k ){
00294 numSubDom = v[k]._numDomain ;
00295 load[numSubDom] += loadAtoms ;
00296 for (unsigned int l = 0 ; l < v.size() ; ++l ){
00297 numSubDom1 = v[l]._numDomain ;
00298 if(v[l]._typeZone != Cut_Zone || numSubDom1 == numSubDom){
00299 (numOfSubDom[numSubDom])[numSubDom1] += 1 ;
00300 }
00301 }
00302 }
00303
00304 }
00305
00306 out << "# " << std::endl
00307 << "# Partition Format version 1.0 " << std::endl
00308 << "# GENERATED BY PartitionBase<TYPEOFATOMS>::writeInFile"<< std::endl
00309 << "#"<< std::endl
00310 << "# Zones of a subdomain and associated zone identifiers :" << std::endl
00311 << "# _____________________________________________ "<< std::endl
00312 << "# | | | | | |"<< std::endl
00313 << "# |Cut Zone|Overlap|Kernel Zone|Overlap|Cut Zone|"<< std::endl
00314 << "# | 2 |Zone 1 | 0 |Zone 1 | 2 |"<< std::endl
00315 << "# |________|_______|___________|_______|________|"<< std::endl
00316 << "# " << std::endl
00317 << "# " << std::endl
00318 << "# " << std::endl
00319 <<"NUMBER OF SUBDOMAINS "<< std::endl << _numberOfPartition << std::endl
00320 <<"RADIUS BUFFER1 " << std::endl << _Rbuffer1 << std::endl
00321 <<"RADIUS BUFFER2 " << std::endl << _Rbuffer2 << std::endl
00322 <<"PARTITION TYPE" << std::endl << _partitionType << std::endl
00323 <<"NUMBER OF ATOMS" << std::endl << _setOfAtoms->numberOfAtoms() << std::endl
00324 << "# " << std::endl ;
00325
00326 out << "# num1 = numero of subdomain" << std::endl
00327 << "# num2 = number of atom(s) in the subdomain num1" << std::endl
00328 << "# num3 = number of AO in the subdomain num1" << std::endl
00329 << "# num4 = number of subdomain(s) overlapping num1" << std::endl
00330 << "# num5 = numero of a subdomain overlapping num1" << std::endl
00331 << "# num6 = number of atom(s) from num1 overlapped by num5"<< std::endl
00332 << "# The doublet (num5, num6) is repeated 'num4-1' time."<< std::endl
00333 << "#num1 num2 num3 num4 num5 num6 " << std::endl
00334 << "PARTITION DESCRIPTION" << std::endl ;
00335 for( nsd = 0 ; nsd < _numberOfPartition ; ++nsd) {
00336
00337
00338 numSubDom = 0;
00339 for ( i = 0 ; i < _numberOfPartition ; ++i){
00340 if( (numOfSubDom[nsd])[i] != 0){
00341 numSubDom++;
00342 }
00343 }
00344 numSubDom--;
00345 out << nsd << " \t"<< _subDomains[nsd].getNumberOfAtoms()
00346 <<" \t" << load[nsd] <<" \t" << numSubDom << " ";
00347
00348 for( i = 0 ; i < _numberOfPartition ; ++i){
00349 if ( (numOfSubDom[nsd])[i] != 0 && i!=nsd){
00350 out << i << " "<< (numOfSubDom[nsd])[i] << " ";
00351 }
00352 }
00353
00354
00355
00356
00357
00358 out << std::endl ;
00359 }
00360 out << "#"<< std::endl
00361 << "# num1 = index of atom"<< std::endl
00362 << "# num2 = the atom num1 is in num2 subdomain(s)"<< std::endl
00363 << "# num3 = numero of a subdomain including num1"<< std::endl
00364 << "# num4 = identifier of zone in num3 where is located num1"<< std::endl
00365 << "# num5 = local numero of num1 in this subdomain num3"<< std::endl
00366 << "# The tripplet (num3, num4, num5) is repeated 'num2-1' time"<< std::endl
00367 << "#"<< std::endl
00368 << "#num1 num2 num3 num4 num5"<< std::endl
00369 << "PARTITION"<< std::endl ;
00370 for(i = 0 ; i < _setOfAtoms->numberOfAtoms() ; ++i) {
00371 out << " "<< i << " "<< _atomsInZone[i].size() ;
00372 std::vector<AtomIn> & v = _atomsInZone[i] ;
00373 for (unsigned int j = 0 ; j < v.size() ; ++j ){
00374 out << " "<< v[j]._numDomain <<" "<< v[j]._typeZone <<" "<< v[j]._localNum ;
00375 }
00376 out << std::endl ;
00377 }
00378 }
00379
00380
00381
00382 template <class TYPEOFATOMS>
00383 void PartitionBase<TYPEOFATOMS>::readFromFile(std::ifstream & data)
00384 {
00385 std::string type,line("") ;
00386 std::stringstream ss ;
00387 int pos, loc =-1,version ;
00388
00389 pos = data.tellg() ;
00390 while ( loc < 0) {
00391 data.seekg(pos);
00392 std::getline(data, line);
00393 loc = line.find_last_of('.');
00394 }
00395 line.erase(line.begin(), line.begin()+loc-1);
00396 ss <<line ; ss >> version ;
00397 data.seekg(pos); ;
00398
00399 if(version == 2) {
00400 std::cout << "newReadFromFile"<<std::endl;
00401 this->newReadFromFile(data);
00402 }
00403 else if(version == 1) {
00404 std::cout << "oldReadFromFile 1"<<std::endl;
00405 this->oldReadFromFile(data);
00406 }
00407 else {
00408 std::cerr << " Wrong format of partition file " << std::endl;
00409 exit(EXIT_FAILURE) ;
00410 }
00411 }
00412 template <class TYPEOFATOMS>
00413 void PartitionBase<TYPEOFATOMS>::oldReadFromFile(std::ifstream & data) {
00414 std::cout << "BEGIN aoldReadFromFile "<<std::endl;
00415 int nbAtoms;
00416 std::string line,comment("#") ;
00417 std::stringstream temp ;
00418
00419 int pos=0 ;
00420 data.seekg(std::ios::beg) ;
00421
00422
00423
00424 std::getline(data,line) ;
00425 while(line[0] == comment[0]) {
00426 std::getline(data,line) ;
00427 }
00428 std::getline(data,line) ;
00429 temp.clear() ;
00430
00431 temp << line ; temp >> _numberOfPartition ; temp.clear();
00432 std::cout << _numberOfPartition <<std::endl;
00433
00434 std::getline(data,line) ;
00435 std::getline(data,line) ;
00436 temp << line ; temp >> _Rbuffer1 ; temp.clear();
00437
00438 std::getline(data,line) ;
00439 std::getline(data,line) ;
00440 temp << line ; temp >> _Rbuffer2 ; temp.clear();
00441
00442 std::getline(data,line) ;
00443 std::getline(data,line) ;
00444 temp << line ; temp >> _partitionType ; temp.clear();
00445
00446 std::getline(data,line) ;
00447 std::getline(data,line) ;
00448 temp << line ; temp >> nbAtoms ; temp.clear();
00449
00450
00451
00452 _subDomains = new subDomain[_numberOfPartition] ;
00453
00454 std::getline(data,line) ;
00455 while(line[0] == comment[0]) {
00456 std::cout << line <<std::endl;
00457 std::getline(data,line) ;
00458 }
00459 std::cout << line <<std::endl;
00460 std::getline(data,line) ;
00461 std::cout << line <<std::endl;
00462
00463
00464
00465 int nSD, nbAtnSD, nbAonSD, nbNeigh, numNeigh, nbATinBuffer1 ;
00466 std::vector< std::vector<int> > voisins(_numberOfPartition) ;
00467 for( int nsd = 0 ; nsd < _numberOfPartition ; ++nsd ){
00468 temp << line ;
00469 temp >> nSD >> nbAtnSD >> nbAonSD >> nbNeigh ;
00470 voisins[nSD].resize(nbNeigh) ;
00471 std::cout << nSD << " "<< nbAtnSD << " "<<nbAonSD << " "<< nbNeigh ;
00472 _subDomains[nsd].initNumberOfOverlapDomain(nbNeigh);
00473 for( int n = 0 ; n < nbNeigh ; ++n ){
00474 temp >> numNeigh >> nbATinBuffer1 ;
00475 voisins[nSD][n] = numNeigh ;
00476 _subDomains[nsd].setNumberOfNeighbor(n, numNeigh) ;
00477 std::cout << " "<<numNeigh << " "<< nbATinBuffer1 ;
00478 }
00479 std::cout << std::endl;
00480 temp.clear();
00481 pos = data.tellg() ;
00482 std::getline(data,line) ;
00483 }
00484 while(line[0] == comment[0]) {
00485 pos = data.tellg() ;
00486 std::getline(data,line) ;
00487 }
00488 data.seekg(pos) ;
00489 std::cout << "Attention lecture des atomes " <<nbAtoms << std::endl ;
00490 int a, aInNBSD, isInSD, typeZone, localNum ;
00491 for( int at = 0 ; at < nbAtoms ; ++at ){
00492 data >> a >> aInNBSD ;
00493 std::cout << a << " "<< aInNBSD << " - ";
00494 for ( int i = 0 ; i < aInNBSD; ++i) {
00495 data >> isInSD>> typeZone >> localNum ;
00496 if( typeZone==0){
00497 _subDomains[isInSD].addIndexInKernel(localNum, at) ;}
00498 else if( typeZone==1){
00499
00500 }
00501 else {
00502
00503 }
00504 std::cout << " "<< isInSD << " "<< typeZone << " "<< localNum ;
00505 }
00506 std::cout << std::endl;
00507 }
00508 std::cerr << "oldFromFile not yet implemented" <<std::endl;
00509
00510
00511 }
00512 template <class TYPEOFATOMS>
00513 void PartitionBase<TYPEOFATOMS>::newReadFromFile(std::ifstream & data)
00514 {
00515 std::stringstream temp ;
00516 int loc ;
00517 std::string line, bid ,comment("#");
00518
00519 std::getline(data, line);
00520 while ( line[0] == comment[0]) {
00521 std::getline(data, line);
00522 }
00523
00524 loc = line.find(":") ;
00525 temp << line.substr(loc+1,line.size() ) ; temp >> _numberOfPartition; temp.clear();
00526 std::getline(data,line) ;
00527 loc = line.find(":") ;
00528 temp << line.substr(loc+1,line.size() ) ; temp >> _Rbuffer1; temp.clear();
00529 std::getline(data,line) ;
00530 loc = line.find(":") ;
00531 temp << line.substr(loc+1,line.size() ) ; temp >> _Rbuffer2; temp.clear();
00532 std::getline(data,line) ;
00533 loc = line.find(":") ;
00534 temp << line.substr(loc+1,line.size() ) ; temp >> _partitionType; temp.clear();
00535
00536
00537 _subDomains = new subDomain[_numberOfPartition] ;
00538
00539 int num ;
00540 for(int i = 0; i < _numberOfPartition ; ++i){
00541
00542 bool cont = true;
00543 do {
00544 std::getline(data,line) ;
00545 loc = line.find("Domain");
00546 if( line.size() !=0 && line.find("Domain") != 2)
00547 { cont = false;}
00548 } while(cont);
00549 loc = line.find(":");
00550 temp << line.substr(loc+1,line.size() ) ; temp >> num; temp.clear();
00551
00552
00553 int nbVoisins, atomInCore, atomInOverlap, atomInCut, numV ;
00554 data >> nbVoisins >> atomInCore>> atomInOverlap >> atomInCut ;
00555 int numAt ;
00556 std::getline(data,line) ;
00557 std::getline(data,line) ;
00558 for(int j = 0; j < atomInCore ; ++j){
00559 data >> numAt ;
00560 _subDomains[num].addIndexInKernel(numAt);
00561 }
00562 std::getline(data,line) ;
00563 std::getline(data,line) ;
00564
00565
00566 _subDomains[num].initNumberOfOverlapDomain(nbVoisins);
00567 for ( int n = 0 ; n < nbVoisins ; ++n){
00568 data >> numV >>atomInOverlap;
00569 _subDomains[num].setNumberOfNeighbor(n, numV) ;
00570 for( int j = 0; j < atomInOverlap ; ++j){
00571 data >> numAt ;
00572 _subDomains[num].addIndexInOverlap(n,numAt);
00573 }
00574 std::getline(data,line) ;
00575 }
00576 std::getline(data,line) ;
00577 for ( int n = 0 ; n < nbVoisins ; ++n){
00578 data >> numV >>atomInCut;
00579 for( int j = 0; j < atomInCut ; ++j){
00580 data >> numAt ;
00581 _subDomains[num].addIndexInCut(n,numAt);
00582 }
00583 std::getline(data,line) ;
00584 }
00585
00586 }
00587
00588 }
00589
00590 template <class TYPEOFATOMS>
00591 void PartitionBase<TYPEOFATOMS>::newWriteInFile(std::ofstream &out)
00592 {
00593
00594
00595
00596 unsigned int num, numOverlap, numCut ;
00597
00598 out <<"# Partition Format Version 2.0 BUG BUG " << std::endl
00599 << "# GENERATED BY PartitionBase<TYPEOFATOMS>::newWriteInFile"<< std::endl ;
00600 out << "NUMBER OF PARTITION : " << _numberOfPartition << std::endl ;
00601 out << "Radius of Buffer1 : " << _Rbuffer1 << std::endl
00602 << "Radius of Buffer2 : " << _Rbuffer2 << std::endl
00603 << "Type : " << _partitionType<< std::endl ;
00604
00605
00606 for(int i = 0; i < _numberOfPartition ; ++i){
00607 numOverlap = numCut = 0 ;
00608 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00609 std::vector<int> & indexOVERLAP = _subDomains[i].getAtomsIndexInOverlap(n) ;
00610 numOverlap += indexOVERLAP.size() ;
00611 std::vector<int> & indexCUT = _subDomains[i].getAtomsIndexInCut(n) ;
00612 numCut += indexCUT.size() ;
00613 }
00614
00615 num = 0 ;
00616 std::vector<int> & indexCORE = _subDomains[i].getAtomsIndexInKernel() ;
00617 out << "Sub Domain number : " << i <<std::endl ;
00618 out << " "<< _subDomains[i].getNumberOfNeighbors() <<" "
00619 << _subDomains[i].getAtomsIndexInKernel().size() <<" "<< numOverlap << " " << numCut<<std::endl ;
00620 out << " CORE ZONE : " <<std::endl <<" ";
00621 if(indexCORE.size() != 0 ){
00622 for(unsigned int j = 0; j < indexCORE.size() ; ++j){
00623 out << indexCORE[j] << " " ;
00624 }
00625 out << std::endl ;
00626 }
00627 num += indexCORE.size() ;
00628 out << " OVERLAP ZONE : " <<std::endl ;
00629 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00630 std::vector<int> & indexOVERLAP = _subDomains[i].getAtomsIndexInOverlap(n) ;
00631 if (indexOVERLAP.size() ==0) {continue ;}
00632 out << " " << _subDomains[i].getNumberOfNeighbor(n) <<" "<< indexOVERLAP.size()<< " ";
00633 for(unsigned int j = 0; j < indexOVERLAP.size() ; ++j){
00634 out << indexOVERLAP[j] << " ";
00635 }
00636 out << std::endl;
00637 num += indexOVERLAP.size() ;
00638 }
00639 out << " CUT ZONE : " <<std::endl ;
00640 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00641 std::vector<int> & indexCUT = _subDomains[i].getAtomsIndexInCut(n) ;
00642 if (indexCUT.size() ==0) {continue ;}
00643 out << " " << _subDomains[i].getNumberOfNeighbor(n) << " "<< indexCUT.size() <<" ";
00644 for(unsigned int j = 0; j < indexCUT.size() ; ++j){
00645 out << indexCUT[j] << " ";
00646 }
00647 out << std::endl;
00648 num += indexCUT.size() ;
00649 }
00650 out << std::endl << std::endl ;
00651 }
00652 }
00653 template <class TYPEOFATOMS>
00654 void PartitionBase<TYPEOFATOMS>::writeInVTKFormat(std::ofstream &out)
00655 {
00656
00657
00658
00659 int nbAtoms = this->_setOfAtoms->numberOfAtoms() ;
00660 std::vector<int> scalar(nbAtoms,-1) ;
00661 for(int i = 0; i < _numberOfPartition ; ++i){
00662 std::vector<int> & indexCORE = _subDomains[i].getAtomsIndexInKernel() ;
00663 for(unsigned int j = 0; j < indexCORE.size() ; ++j){
00664 scalar[indexCORE[j]] = i ;
00665 }
00666 if(_partitionType == "QC_Dixon_Part"){continue;}
00667
00668 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00669 std::vector<int> & indexOVERLAP = _subDomains[i].getAtomsIndexInOverlap(n) ;
00670 if( i > _subDomains[i].getNumberOfNeighbor(n) ){
00671 for(unsigned int j = 0; j < indexOVERLAP.size() ; ++j){
00672 scalar[indexOVERLAP[j]] = _numberOfPartition ;
00673 }
00674 }
00675 }
00676 }
00677
00678 out << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<< std::endl
00679 <<"<PolyData>" << std::endl <<"<Piece NumberOfPoints=\"" << nbAtoms <<"\">"<< std::endl ;
00680 out << "<Points>" <<std::endl ;
00681 out << "<DataArray type=\"Float64\" Name=\"coord_atoms\" NumberOfComponents=\"3\" format=\"ascii\">" <<std::endl ;
00682 double * xyz ;
00683 for ( int i = 0 ; i < nbAtoms ; ++i){
00684 xyz = this->_setOfAtoms->getElement(i).pos();
00685 out << xyz[0] << " " <<xyz[1] << " " <<xyz[2] << " " <<std::endl;
00686 }
00687 out << "</DataArray>" <<std::endl ;
00688 out << "</Points>" <<std::endl ;
00689
00690 out << "<PointData Scalars=\"PartitionNumber\"> "<<std::endl ;
00691 out << "<DataArray type=\"Int32\" Name=\"PartitionNumber\" format=\"ascii\"> " <<std::endl;
00692
00693 for ( unsigned int i = 0 ; i < scalar.size() ; ++i){
00694 out << scalar[i] << " " ;
00695 }
00696 out <<std::endl ;
00697
00698 out << "</DataArray>" <<std::endl ;
00699 out << "</PointData>" <<std::endl ;
00700 out << "</Piece>" <<std::endl ;
00701 out << "</PolyData>" <<std::endl ;
00702 out << "</VTKFile>" <<std::endl ;
00703
00704 }
00705
00706
00707 template <class TYPEOFATOMS>
00708 void PartitionBase<TYPEOFATOMS>::newWriteInFileOld(std::ofstream &out)
00709 {
00710
00711
00712
00713 unsigned int num, numOverlap, numCut ;
00714
00715 out << "NEW FORMAT" << std::endl ;
00716 out << "NUMBER OF PARTITION : " << _numberOfPartition << std::endl ;
00717 out << "Size of Buffer1 : " << _Rbuffer1 << std::endl
00718 << "Size of Buffer2 : " << _Rbuffer2 << std::endl
00719 << "Type : " << _partitionType<< std::endl ;
00720
00721
00722 for(int i = 0; i < _numberOfPartition ; ++i){
00723 numOverlap = numCut = 0 ;
00724 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00725 std::vector<int> & indexOVERLAP = _subDomains[i].getAtomsIndexInOverlap(n) ;
00726 numOverlap += indexOVERLAP.size() ;
00727 std::vector<int> & indexCUT = _subDomains[i].getAtomsIndexInCut(n) ;
00728 numCut += indexCUT.size() ;
00729 }
00730
00731 num = 0 ;
00732 std::vector<int> & indexCORE = _subDomains[i].getAtomsIndexInKernel() ;
00733 out << "Sub Domain number : " << i <<std::endl ;
00734 out << " "<< _subDomains[i].getAtomsIndexInKernel().size() <<" "<< numOverlap << " " << numCut<<std::endl ;
00735 out << " CORE ZONE : " ;
00736 for(unsigned int j = 0; j < indexCORE.size() ; ++j){
00737 out << " { " << j << " , "<< indexCORE[j] << " } ";
00738 }
00739 out << std::endl ;
00740 num += indexCORE.size() ;
00741 out << " OVERLAP ZONE : " ;
00742 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00743 std::vector<int> & indexOVERLAP = _subDomains[i].getAtomsIndexInOverlap(n) ;
00744 for(unsigned int j = 0; j < indexOVERLAP.size() ; ++j){
00745 out << " { " << num + j << " , "<< indexOVERLAP[j] << " } ";
00746 }
00747 num += indexOVERLAP.size() ;
00748 }
00749 out << std::endl ;
00750 out << " CUT ZONE : " ;
00751 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00752 std::vector<int> & indexCUT = _subDomains[i].getAtomsIndexInCut(n) ;
00753 for(unsigned int j = 0; j < indexCUT.size() ; ++j){
00754 out << " { " << num + j << " , "<< indexCUT[j] << " } ";
00755 }
00756 num += indexCUT.size() ;
00757 }
00758 out << std::endl << std::endl ;
00759 out <<"======================================================"<<std::endl
00760 << std::endl << std::endl ;
00761 }
00762 }
00763 template <class TYPEOFATOMS>
00764 void PartitionBase<TYPEOFATOMS>::getAtomsIndexInPartition(const int& part, std::vector<int>& index){
00765 std::cout << _subDomains[part].getNumberOfAtoms() <<std::endl;
00766 index.resize(_subDomains[part].getNumberOfAtoms() );
00767 std::vector<int> & indexCORE = _subDomains[part].getAtomsIndexInKernel() ;
00768 for(unsigned int j = 0; j < indexCORE.size() ; ++j){
00769 index[j] = indexCORE[j] ;
00770 }
00771 int num = indexCORE.size();
00772 for ( int n = 0 ; n < _subDomains[part].getNumberOfNeighbors() ; ++n){
00773 std::vector<int> & indexOVERLAP = _subDomains[part].getAtomsIndexInOverlap(n) ;
00774 for(unsigned int j = 0; j < indexOVERLAP.size() ; ++j){
00775 index[num + j] = indexOVERLAP[j] ;
00776 }
00777 num += indexOVERLAP.size() ;
00778 }
00779 for ( int n = 0 ; n < _subDomains[part].getNumberOfNeighbors() ; ++n){
00780 std::vector<int> & indexCUT = _subDomains[part].getAtomsIndexInCut(n) ;
00781 for(unsigned int j = 0; j < indexCUT.size() ; ++j){
00782 index[num + j] = indexCUT[j] ;
00783 }
00784 num += indexCUT.size() ;
00785 }
00786 }
00787
00788
00789
00790
00791 template <class TYPEOFATOMS>
00792 void PartitionBase<TYPEOFATOMS>::buildAtomsInZone() {
00793
00794
00795 AtomIn atomNSD ;
00796 int num ;
00797
00798 for (unsigned int i = 0 ; i < _atomsInZone->size() ; ++i){
00799 _atomsInZone[i].clear() ;
00800 }
00801
00802
00803 for(int i = 0; i < _numberOfPartition ; ++i){
00804
00805
00806 num = 0 ;
00807 atomNSD._numDomain = i ;
00808
00809
00810
00811 std::vector<int> & indexCORE = _subDomains[i].getAtomsIndexInKernel() ;
00812 atomNSD._typeZone = Core_Zone ;
00813 for(unsigned int j = 0; j < indexCORE.size() ; ++j){
00814 atomNSD._localNum = j ;
00815 _atomsInZone[indexCORE[j]].push_back(atomNSD);
00816 }
00817
00818
00819
00820 num += indexCORE.size() ;
00821 atomNSD._typeZone = Overlap_Zone ;
00822 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00823
00824 std::vector<int> & indexOVERLAP = _subDomains[i].getAtomsIndexInOverlap(n) ;
00825 for(unsigned int j = 0; j < indexOVERLAP.size() ; ++j){
00826 atomNSD._localNum = num + j ;
00827 _atomsInZone[indexOVERLAP[j]].push_back(atomNSD);
00828 }
00829 num += indexOVERLAP.size() ;
00830 }
00831
00832
00833
00834 atomNSD._typeZone = Cut_Zone ;
00835 for ( int n = 0 ; n < _subDomains[i].getNumberOfNeighbors() ; ++n){
00836
00837
00838 std::vector<int> & indexCUT = _subDomains[i].getAtomsIndexInCut(n) ;
00839 for(unsigned int j = 0; j < indexCUT.size() ; ++j){
00840 atomNSD._localNum = num + j ;
00841
00842 _atomsInZone[indexCUT[j]].push_back(atomNSD);
00843 }
00844
00845 num += indexCUT.size() ;
00846 }
00847
00848 }
00849 }
00850 template <class TYPEOFATOMS>
00851 bool PartitionBase<TYPEOFATOMS>::setAtomsInOverlapAndCut(){
00852
00853 bool retVal = false;
00854
00855 int k, nsd, nsd_k = 0, atomJ, fragNum ;
00856 double dist=0., distMin ;
00857 double R = _Rbuffer1 + _Rbuffer2, R1 = _Rbuffer1* _Rbuffer1 ;
00858
00859 std::vector<bool> alreadyUsed ;
00860
00861 R = R*R ;
00862
00863 std::vector<int> numFragInSD,numFragInSD_k ;
00864
00865 for(nsd = _numberOfPartition -1 ; nsd > -1; --nsd) {
00866
00867
00868 std::vector<int>& index = _subDomains[nsd].getAtomsIndexInKernel();
00869
00870
00871
00872 bool isIn ;
00873 if( _useFrag){
00874 numFragInSD.clear() ;
00875 for ( unsigned int i = 0 ; i < index.size() ; ++i){
00876 isIn = false ;
00877 fragNum = _setOfAtoms->getElement(index[i]).getFragmentNumber() -1;
00878 for ( unsigned int k = 0 ; k < numFragInSD.size() ; ++k){
00879 if(fragNum == numFragInSD[k]){ isIn = true ; break ;}
00880 }
00881 if(!isIn) {
00882 numFragInSD.push_back(fragNum) ;
00883 }
00884 }
00885 }
00886
00887 for (k = 0 ; k < _subDomains[nsd].getNumberOfNeighbors() ; ++k){
00888 nsd_k = _subDomains[nsd].getNumberOfNeighbor(k) ;
00889
00890 std::vector<int>& indexZONE1 = _subDomains[nsd].getAtomsIndexInOverlap(k);
00891 std::vector<int>& indexZONE2 = _subDomains[nsd].getAtomsIndexInCut(k);
00892 indexZONE1.clear();
00893 indexZONE2.clear();
00894
00895 std::vector<int>& index_k = _subDomains[nsd_k].getAtomsIndexInKernel();
00896
00897
00898
00899
00900
00901 if( !_useFrag){
00902 for ( unsigned int j = 0 ; j < index_k.size() ; ++j){
00903 atomJ = index_k[j] ;
00904 distMin = 9999. ;
00905 for ( unsigned int i = 0 ; i < index.size() ; ++i){
00906 dist = this->_setOfAtoms->getElement(atomJ).dist2To(this->_setOfAtoms->getElement(index[i]).pos() );
00907 if (distMin > dist){
00908 distMin = dist ;
00909 }
00910 }
00911 if( distMin > R){
00912 continue ; }
00913 else {
00914 if (distMin < R1 ){
00915 indexZONE1.push_back(atomJ) ;
00916 }
00917 else {
00918 indexZONE2.push_back( atomJ) ;
00919 }
00920 }
00921 }
00922 }
00923 else {
00924 std::cout << "Construction avec des fragments " <<std::endl;
00925 alreadyUsed.resize(index_k.size()) ;
00926 numFragInSD_k.clear() ;
00927 for ( unsigned int j = 0 ; j < index_k.size() ; ++j){
00928 alreadyUsed[j] = false ;
00929 isIn = false ;
00930 fragNum = _setOfAtoms->getElement(index_k[j]).getFragmentNumber() -1;
00931 for ( unsigned int k = 0 ; k < numFragInSD_k.size() ; ++k){
00932 if(fragNum == numFragInSD_k[k]){ isIn = true ; break ;}
00933 }
00934 if(!isIn) {
00935 numFragInSD_k.push_back(fragNum) ;
00936 }
00937 }
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951 int atIdx, atIdx_k, fragNum,fragNum_k,atIdxToAdd ;
00952 bool isInShell1, isInShell2 ;
00953 for ( unsigned int frag = 0 ; frag < numFragInSD.size() ; ++frag){
00954 fragNum = numFragInSD[frag] ;
00955 atIdx = _firstATofFrag[fragNum] ;
00956
00957 for ( unsigned int frag_k = 0 ; frag_k < numFragInSD_k.size() ; ++frag_k){
00958 isInShell1 = isInShell2 = false ;
00959 fragNum_k = numFragInSD_k[frag_k] ;
00960
00961 while ( atIdx != -1 && !isInShell1){
00962 atIdx_k = _firstATofFrag[fragNum_k] ;
00963
00964
00965 while ( atIdx_k != -1 && !isInShell1){
00966 dist = this->_setOfAtoms->getElement(atIdx).dist2To(this->_setOfAtoms->getElement(atIdx_k).pos() );
00967 if( dist < R){
00968 if( dist < R1){
00969 isInShell1 = true ;
00970 }
00971 else {
00972 isInShell2 = true ;
00973 }
00974 }
00975
00976
00977
00978 atIdx_k = _nextATofFrag[atIdx_k] ;
00979 }
00980
00981 atIdx = _nextATofFrag[atIdx] ;
00982 }
00983 atIdx = _firstATofFrag[fragNum] ;
00984 if(isInShell1){
00985 atIdxToAdd = _firstATofFrag[fragNum_k] ;
00986 while ( atIdxToAdd != -1){
00987 indexZONE1.push_back(atIdxToAdd) ;
00988 atIdxToAdd = _nextATofFrag[atIdxToAdd] ;
00989 }
00990 std::cout << "On ajoute les atomes du fragment "<< frag_k<< " dans le shell 1 du fragment "<<frag <<std::endl;
00991 break ;
00992 }
00993 if(isInShell2){
00994 atIdxToAdd = _firstATofFrag[fragNum_k] ;
00995 while ( atIdxToAdd != -1){
00996 indexZONE2.push_back(atIdxToAdd) ;
00997 atIdxToAdd = _nextATofFrag[atIdxToAdd] ;
00998 }
00999 std::cout << "On ajoute les atomes du fragment "<< frag_k<< " dans le shell 2 du fragment "<<frag <<std::endl;
01000 break ;
01001 }
01002 }
01003 }
01004 }
01005 }
01006 }
01007
01008 return retVal;
01009 }
01010
01011
01012
01013 template class PartitionBase<PDBstructure> ;
01014 template class PartitionBase<setOfAtoms> ;
01015
01016