00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
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
00114
00115
00116
00117
00118
00119
00120
00121 }
00122
00123 template <class TYPEOFSYSTEM>
00124 bool subDomain::addAtomsFromSubdomainInOverlapAndCut(const int & localOverlapSubDomain,
00125 subDomain & domain,
00126 const TYPEOFSYSTEM &atoms,
00127
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
00144 }
00145
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
00166
00167
00168
00169 for(int at = 0 ; at < static_cast<int>(index2.size()) ; ++at){
00170 numAtom = index2[at] ;
00171
00172
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
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
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 std::vector<int>::iterator startIterator = index2.begin();
00206 index2.erase(startIterator,startIterator + indexToRemove.size());
00207
00208
00209
00210
00211
00212
00213 indexToRemove.clear() ;
00214
00215
00216
00217
00218 for(int at = 0 ; at < static_cast<int>(index1.size()) ; ++at){
00219 numAtom = index1[at] ;
00220
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
00226
00227
00228
00229 if ( dist > R )
00230 { continue ;}
00231
00232 if ( dist < R_halfOverlap ){
00233 this->addIndexInOverlap(localOverlapSubDomain1,numAtom) ;
00234 _atomIsInKernel[at] = 0;
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
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 startIterator = index1.end() - indexToRemove.size();
00260 index1.erase(startIterator, index1.end());
00261
00262
00263
00264
00265
00266
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
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
00313 }
00314
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
00326
00327
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
00335 while ( ( localNum < _numberOfAtoms) ) {
00336
00337
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
00347 if (xInOver[imin] > x ){
00348
00349
00350 atomK._localNum = localNum ;
00351 setLocalNum(atomsInZone[_indexInKernel[atK]],atomK) ;
00352 localNum++; atK++;
00353 x = 999999999.0;
00354
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
00366
00367 atomO._localNum = localNum ;
00368 setLocalNum(atomsInZone[(_indexInOverlap[imin])[atOIn[imin]]],atomO );
00369 localNum++;
00370
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
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
00391 setLocalNum(atomsInZone[(_indexInCut[i])[j]],atomO) ;
00392 }
00393 }
00394
00395 }
00396
00397
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) ;