00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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
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 }
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
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
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
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
00157 }
00158
00159
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
00180
00181
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
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
00213
00214 this->innerReperOfMolecule() ;
00215 this->projectOnInnerReper() ;
00216
00217
00218
00219 int* ordre ; ordre = new int [numberOfAtoms] ;
00220 for(i = 0 ; i < numberOfAtoms ; ++i) {
00221 ordre[i] = i;
00222 }
00223 this->sortInFirstDimension(ordre) ;
00224
00225
00226
00227
00228 this->setNeigboors();
00229
00230 if(_useFrag){
00231 build_fragmentStructure( ) ;
00232 }
00233
00234
00235
00236
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
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
00259
00260
00261 }
00262
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] ;
00289
00290
00291
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
00305
00306
00307 int equilibrateCpt;
00308 int nbMaxOfIterationToEquilibrate = numberOfIterations;
00309 for(equilibrateCpt = 0 ; equilibrateCpt < nbMaxOfIterationToEquilibrate ; ++equilibrateCpt) {
00310
00311
00312
00313
00314 this->setAtomsInKernel() ;
00315
00316
00317
00318
00319
00320
00321 this->setAtomsInOverlapAndCut();
00322
00323
00324 for (int nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00325 _subDomains[nsd].computeNumberOfAtoms();
00326 }
00327
00328 this->setLoad() ;
00329
00330
00331 std::cout <<" Iteration : " << equilibrateCpt <<std::endl;
00332 isCV = computeBoxEquilibrateCharge(R_overlap, R_cut);
00333 if(isCV || nbMaxOfIterationToEquilibrate == equilibrateCpt + 1)
00334 {break ;}
00335
00336
00337 this->resetData();
00338 }
00339
00340
00341 this->buildLocalNumerotation();
00342
00343
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
00358
00359
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
00365 _subDomains[nsd].setMyNumber(nsd,nsdx ,nsdy, nsdz);
00366 }
00367 }
00368 }
00369
00370
00371
00372 int kz, kz_deb, kz_fin, ky, ky_deb, ky_fin, kx, kx_deb, kx_fin ;
00373
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
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
00432
00433 for ( nsd = 0 ; nsd < _numberOfPartition ; ++nsd){
00434 _subDomains[nsd].initNumberOfOverlapDomain() ;
00435 }
00436
00437
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
00444
00445
00446
00447
00448
00449
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
00462
00463 j = 0 ;
00464 for ( nsd = 0 ; nsd < _numberOfPartition; ++nsd){
00465
00466
00467
00468
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
00475
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
00487 _subDomains[nsd].addIndexInKernel(j) ;
00488 }
00489
00490 }
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
00508 }
00509 if(numberOfAtoms-nbAtomsInSubDomains!=0)
00510 {std::cerr << "NumberOfAtoms : "<< numberOfAtoms<< " Perte d'atomes : "<< numberOfAtoms-nbAtomsInSubDomains<<std::endl;}
00511
00512
00513
00514 if(emptyKernel.size() !=0 )
00515 {
00516 for ( nsd = 0 ; nsd < _numberOfPartition; ++nsd){
00517
00518 }
00519 }
00520 return retVal ;
00521 }
00522
00523 bool BoxPartition::setAtomsInOverlapAndCut(){
00524 bool retVal = false;
00525
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
00533 for (k = 0 ; k < _subDomains[nsd].getNumberOfNeighbors() ; ++k){
00534 nsd_k = _subDomains[nsd].getNumberOfNeighbor(k) ;
00535 if( nsd_k < nsd )
00536 { continue ; }
00537
00538 _subDomains[nsd].addAtomsFromSubdomainInOverlapAndCut(k, _subDomains[nsd_k],*_setOfAtoms,
00539 R_overlap, R_cut,_useFrag);
00540 }
00541
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
00550 }
00551 }
00552
00553
00554 void BoxPartition::sortInFirstDimension(int* list) {
00555
00556
00557
00558
00559
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 }
00580
00581 }
00582 }
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 }
00609 meanCharge /= _numberOfPartition ;
00610
00611
00612
00613 totalTmpWideSD = 0; nb = 0 ; removeLength = 0.0 ;
00614 for(nsd = 0 ; nsd < _numberOfPartition ; ++nsd) {
00615 tmp = c[nsd]*RhalfOverlapAndCut ;
00616
00617
00618
00619 tmpWideSD[nsd] = (meanCharge / static_cast<double>(_subDomains[nsd].load()))*(wideSD[nsd] + tmp) - tmp;
00620
00621
00622
00623
00624 if(tmpWideSD[nsd] < c[nsd]*RhalfOverlap){
00625 removeLength += c[nsd]*RhalfOverlap ; }
00626 else{
00627 ++nb ; totalTmpWideSD += tmpWideSD[nsd] ;
00628 }
00629 }
00630
00631
00632 gamma = (totalWideSD - removeLength ) / totalTmpWideSD;
00633 totalWideSD = 0.0 ; diff = 0.0; errorLoad = 0 ;
00634
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
00642 errorLoad = std::max( errorLoad , meanCharge - _subDomains[nsd].load()) ;
00643
00644 diff = std::max( diff , fabs(wideSD[nsd] - tmpWideSD[nsd]) );
00645 wideSD[nsd] = tmpWideSD[nsd];
00646 totalWideSD += wideSD[nsd] ;
00647 }
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
00656
00657 if( error < 1.0e-2 )
00658 {isCV = true ;
00659 std::cout << " La methode a convergee." << std::endl;
00660 }
00661
00662
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
00674 }
00675
00676
00677 std::cout << "======================================================================="
00678 << std::endl<< std::endl;
00679
00680 delete[] wideSD ; wideSD = NULL;
00681 delete[] tmpWideSD; tmpWideSD = NULL;
00682 return isCV;
00683 }