00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <fstream>
00022 #include <vector>
00023
00024
00025 #include "QCCommon.hpp"
00026 #include "QCTable.hpp"
00027 #include "QCPartitioner.hpp"
00028 #include "QCGeneralData.hpp"
00029 #include "QCSystem.hpp"
00030 #include "QCGlobalSystem.hpp"
00031 #include "QCMDSystem.hpp"
00032 #include "QCSCFAlgo.hpp"
00033 #include "QCDCAlgo.hpp"
00034
00035
00036
00037
00041 template <class TPSystem>
00042 QCPartitioner<TPSystem>::QCPartitioner (void) : isReady(false), readPartFile(true), useFrag(false),
00043 nbPartitions(0), partType(QC_STANDARD_PART), radius1(3.0),
00044 radius2(4.0), charge(NULL), _ordering(NULL), files(NULL)
00045 {}
00046
00050 template <class TPSystem>
00051 QCPartitioner<TPSystem>::~QCPartitioner (void) {
00052 QC_TRACE_END("BEGIN QCPartitioner<TPSystem>::~QCPartitioner");
00053 if (this->charge != NULL) {
00054 delete [] this->charge; this->charge = NULL;
00055 }
00056 this->isReady = false ;
00057 QC_TRACE_END("END QCPartitioner<TPSystem>::~QCPartitioner");
00058 }
00059
00060
00061 template <class TPSystem>
00062 void
00063 QCPartitioner<TPSystem>::initData( const QCGeneralData& data, const QCFiles& files_){
00064 QC_TRACE_INIT("BEGIN QCPartitioner<TPSystem>::initData");
00065
00066 this->readPartFile = data.isReadPartitionFromFile();
00067 this->nbPartitions = data.getNbPartitions();
00068 this->partType = data.getPartitionType();
00069 this->partTypeStr = data.getPartitionTypeStr();
00070 this->radius1 = data.getBuffDiameter();
00071 this->radius2 = data.getCutRadius();
00072 this->isReady = true ;
00073 this->files = &files_;
00074
00075 QC_TRACE_INIT("END QCPartitioner<TPSystem>::initData");
00076 }
00080 template <class TPSystem>
00081 void
00082 QCPartitioner<TPSystem>::init (const QCGeneralData& data, const QCFiles& files_,
00083 const int * load, const int sizeLoad) {
00084
00085 QC_TRACE_INIT("BEGIN QCPartitioner<TPSystem>::init");
00086
00087
00088 readPartFile = data.isReadPartitionFromFile();
00089 if(!this->readPartFile){
00090 this->nbPartitions = data.getNbPartitions();
00091 this->partType = data.getPartitionType();
00092 this->partTypeStr = data.getPartitionTypeStr();
00093 this->radius1 = data.getBuffDiameter();
00094 this->radius2 = data.getCutRadius();
00095 this->isReady = true ;
00096 }
00097 else {
00098 this->isReady = false ;
00099 this->partTypeStr = "NONE";
00100 }
00101 this->files = &files_;
00102 this->setLoad(load,sizeLoad);
00103 QC_TRACE_INIT("END QCPartitioner<TPSystem>::init");
00104 }
00105
00106 template <class TPSystem>
00107 void QCPartitioner<TPSystem>::setLoad(const int * load, const int sizeLoad){
00108
00109 QC_TRACE_INIT("BEGIN QCPartitioner<TPSystem>::setLoad(" <<load<<", "<<sizeLoad<<")");
00110 if (load != NULL && sizeLoad != 0) {
00111 if(this->charge != NULL)
00112 { delete [] this->charge;}
00113 QC_TRACE_INIT(" QCPartitioner<TPSystem>::setLoad(" << this->charge<<")");
00114 this->charge = new int [sizeLoad];
00115
00116 for (int i=0; i < sizeLoad; ++i) {
00117 this->charge[i] = load[i];
00118
00119 }
00120 }
00121 QC_TRACE_INIT("END QCPartitioner<TPSystem>::setLoad");
00122 }
00126 template <class TPSystem>
00127 void
00128 QCPartitioner<TPSystem>::partitioning (TPSystem& ) {
00129
00130 cout << "virtual QCPartitioner::partitioning(): NOTHING TO DO HERE." << endl;
00131 exit(EXIT_FAILURE) ;
00132 }
00136 template <class TPSystem>
00137 void
00138 QCPartitioner<TPSystem>::readFromFile (TPSystem& system, QCGeneralData& data,
00139 const QCFiles& Files,const string& path) {
00140
00141 QC_TRACE_INIT("BEGIN QCPartitioner<TPSystem>::readFromFile");
00142
00143 QCReader reader;
00144 std::string name ;
00145
00146 name = Files.getPartitionFile() ;
00147
00148 reader.readPartitionFile(name, data, system, path);
00149
00150 this->initData(data,Files);
00151
00152 for( int nsd = 0 ; nsd < this->nbPartitions ; ++nsd) {
00153
00154 system.getSubDomain(nsd).updateInternalStructure() ;
00155 }
00156
00157
00158
00159 QC_TRACE_INIT("END QCPartitioner<TPSystem>::readFromFile");
00160 }
00165 template <class TPSystem>
00166 void QCPartitioner<TPSystem>::updateInternalStructure(TPSystem& system){
00168 QC_TRACE_INIT("BEGIN QCPartitioner<TPSystem>::updateInternalStructure "<< this->nbPartitions) ;
00169
00170 int nsd ;
00171 QC_TRACE_INIT(" nbPartitions "<<this->nbPartitions<<" system.getNbDomains() : "<< system.getNbDomains());
00172
00173 for( int nsd = 0 ; nsd < system.getNbDomains() ; ++nsd) {
00174
00175 system.getSubDomain(nsd).updateInternalStructure() ;
00176 }
00177 QC_TRACE_INIT("END QCPartitioner<TPSystem>::updateInternalStructure") ;
00178
00179 }
00180
00181 template <class TPSystem>
00182 void QCPartitioner<TPSystem>::writeInFile (const TPSystem& system, const std::string& fileNameIn,
00183 const string& path)
00184 {
00185 QC_TRACE_OUT("BEGIN QCPartitioner<TPSystem>::writeInFile");
00186 std::ofstream out;
00187 std::string fileName;
00188
00189 fileName = path +"/"+ fileNameIn + ".part1";
00190 QC_TRACE_OUT(" write partition in Old QCPP format in file : " << fileName );
00191 out.open(fileName.c_str());
00192 this->writeInFileOLD(out,system);
00193 out.close();
00194
00195 fileName = path +"/"+ fileNameIn + ".part2";
00196 out.open(fileName.c_str());
00197 QC_TRACE_OUT(" write partition in New QCPP format in file : " << fileName );
00198 this->writeInFileNew(out,system);
00199 out.close();
00200
00201 fileName = fileNameIn + "-part.vtp";
00202
00203 QC_TRACE_OUT(" write partition in VTK format in file : " << fileName );
00204 out.open(fileName.c_str());
00205 this->writeInFileVTK(out, system );
00206
00207 QC_TRACE_OUT("END QCPartitioner<TPSystem>::writeInFile");
00208 }
00209 template <class TPSystem>
00210 void QCPartitioner<TPSystem>::writeHeader ( std::ofstream & out)
00211 {
00212 QC_TRACE_OUT("BEGIN QCPartitioner<TPSystem>::writeHeader");
00213
00214 out << "NUMBER OF PARTITION : " << this->nbPartitions << std::endl ;
00215 out << "Radius of Buffer1 : " << this->radius1 << std::endl
00216 << "Radius of Buffer2 : " << this->radius2 << std::endl
00217 << "Type : " << this->partTypeStr<< std::endl ;
00218 QC_TRACE_OUT("END QCPartitioner<TPSystem>::writeHeader");
00219 }
00220
00221
00222 template <class TPSystem>
00223 void QCPartitioner<TPSystem>::writeInFileNew ( std::ofstream & out,const TPSystem& system)
00224 {
00225 QC_TRACE_OUT("BEGIN QCPartitioner<TPSystem>::writeInFileNew");
00226
00227
00228 out << "# Partition Format version 2.0" << std::endl
00229 << "# GENERATED BY QCPartitioner<TPSystem>::writeInFileV2"<< std::endl;
00230 this->writeHeader(out) ;
00231
00232 for(int i = 0; i < this->nbPartitions ; ++i){
00233 system.getSubDomain(i).writeInFileNew(out) ;
00234 }
00235 out << "END "<<std::endl;
00236 QC_TRACE_OUT("END QCPartitioner<TPSystem>::writeInFileNew");
00237 }
00238
00239 template <class TPSystem>
00240 void QCPartitioner<TPSystem>::writeInFileOLD ( std::ofstream & out,const TPSystem& system)
00241 {
00242 QC_TRACE_OUT("BEGIN QCPartitioner<TPSystem>::writeInFileOLD");
00243
00244 out << "# Partition Format Version 1.0" << std::endl
00245 << "# GENERATED BY QCPartitioner<TPSystem>::writeInFileV1"<< std::endl
00246 << "#"<< std::endl
00247 << "# Zones of a subdomain and associated zone identifiers :" << std::endl
00248 << "# _____________________________________________ "<< std::endl
00249 << "# | | | | | |"<< std::endl
00250 << "# |Cut Zone|Overlap|Kernel Zone|Overlap|Cut Zone|"<< std::endl
00251 << "# | 2 |Zone 1 | 0 |Zone 1 | 2 |"<< std::endl
00252 << "# |________|_______|___________|_______|________|"<< std::endl
00253 << "# " << std::endl
00254 << "# " << std::endl
00255 << "# " << std::endl
00256
00257 <<"NUMBER OF SUBDOMAINS "<< std::endl << this->nbPartitions << std::endl
00258 <<"RADIUS BUFFER1 " << std::endl << this->radius1 << std::endl
00259 <<"RADIUS BUFFER2 " << std::endl << this->radius2 << std::endl
00260 <<"PARTITION TYPE" << std::endl << this->partTypeStr << std::endl
00261 <<"NUMBER OF ATOMS" << std::endl << system.getNbAtoms() << std::endl
00262 << "# " << std::endl ;
00263
00264 out << "# num1 = numero of subdomain" << std::endl
00265 << "# num2 = number of atom(s) in the subdomain num1" << std::endl
00266 << "# num3 = number of AO in the subdomain num1" << std::endl
00267 << "# num4 = number of subdomain(s) overlapping num1" << std::endl
00268 << "# num5 = numero of a subdomain overlapping num1" << std::endl
00269 << "# num6 = number of atom(s) from num1 overlapped by num5"<< std::endl
00270 << "# The doublet (num5, num6) is repeated 'num4-1' time."<< std::endl
00271 << "#num1 num2 num3 num4 num5 num6 "<< std::endl
00272 << "PARTITION DESCRIPTION" << std::endl ;
00273 int i,nsd ;
00274 std::vector<std::vector<QCAtomIn> > atomsInZone(system.getNbAtoms()) ;
00275 for( nsd = 0 ; nsd < nbPartitions ; ++nsd) {
00276 #ifdef TRACE_OUT
00277 system.getSubDomain(nsd).printHeadOfInternalStructure();
00278 #endif
00279 system.getSubDomain(nsd).writeSubDomainHeader(out) ;
00280 system.getSubDomain(nsd).fillAtomInStructure(atomsInZone);
00281 #ifdef TRACE_OUT
00282 system.getSubDomain(nsd).printHeadOfInternalStructure();
00283 #endif
00284 }
00285
00286 out << "#" <<std::endl;
00287 out << "# num1 = index of atom" << std::endl ;
00288 out << "# num2 = the atom num1 is in num2 subdomain(s)" << std::endl ;
00289 out << "# num3 = numero of a subdomain including num1" << std::endl ;
00290 out << "# num4 = identifier of zone in num3 where is located num1" << std::endl ;
00291 out << "# num5 = local numero of num1 in this subdomain num3" << std::endl ;
00292 out << "# The tripplet (num3, num4, num5) is repeated 'num2-1' time" << std::endl ;
00293 out << "#" << std::endl ;
00294 out << "#num1 num2 num3 num4 num5 " << std::endl ;
00295 out << "PARTITION" << std::endl ;
00296
00297 for( i = 0 ; i < system.getNbAtoms() ; ++i)
00298 {
00299 out << " "<< i << " "<< atomsInZone[i].size() ;
00300 std::vector<QCAtomIn> & v = atomsInZone[i] ;
00301 for (unsigned int j = 0 ; j < v.size() ; ++j ){
00302 out << " "<< v[j].numDomain <<" "<< v[j].typeZone <<" "<< v[j].localNum ;
00303 }
00304 out << std::endl ;
00305 }
00306 QC_TRACE_OUT("END QCPartitioner<TPSystem>::writeInFileOLD");
00307 }
00308
00309 template <class TPSystem>
00310 void QCPartitioner<TPSystem>::writeInFileVTK(std::ofstream & out, const TPSystem& system){
00311 QC_TRACE_OUT("BEGIN QCPartitioner<TPSystem>::writeInFileVTK");
00312
00313 const QCSystem& globalSystem = system.getRootSystem();
00314 int nbAtoms = globalSystem.getNbAtoms();
00315
00316 std::vector<int> scalar(nbAtoms,-1) ;
00317
00318 for(int i = 0; i < this->nbPartitions ; ++i){
00319
00320 QCSubDomain & domaine = system.getSubDomain(i) ;
00321 QCSubDomainZone * zone = system.getSubDomain(i).getZone();
00322 const int * globalIndex = domaine.getGlobalIdxs() ;
00323
00324
00325 for( int j = 0; j < domaine.getNbAtoms() ; ++j){
00326
00327
00328 if(zone[j] == QC_CORE){
00329 scalar[globalIndex[j]] = i ;
00330 }
00331 else if (this->partType == QC_DIXON_PART ) {
00332
00333 continue;
00334 }
00335 else if(zone[j] == QC_SHELL1 ) {
00336 scalar[globalIndex[j]] = this->nbPartitions ;
00337 }
00338
00339 }
00340 }
00341
00342 out << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<< std::endl
00343 <<"<PolyData>" << std::endl <<"<Piece NumberOfPoints=\"" << nbAtoms <<"\">"<< std::endl ;
00344 out << "<Points>" <<std::endl ;
00345 out << "<DataArray type=\"Float64\" Name=\"coord_atoms\" NumberOfComponents=\"3\" format=\"ascii\">" <<std::endl ;
00346 QCPoint3D xyz ;
00347 for ( int i = 0 ; i < nbAtoms ; ++i){
00348 globalSystem.getPointAt(i,xyz);
00349 out << xyz[0] << " " <<xyz[1] << " " <<xyz[2] << " " <<std::endl;
00350 }
00351 out << "</DataArray>" <<std::endl ;
00352 out << "</Points>" <<std::endl ;
00353
00354 out << "<PointData Scalars=\"PartitionNumber\"> "<<std::endl ;
00355 out << "<DataArray type=\"Int32\" Name=\"PartitionNumber\" format=\"ascii\"> " <<std::endl;
00356
00357 for ( unsigned int i = 0 ; i < scalar.size() ; ++i){
00358 out << scalar[i] << " " ;
00359 }
00360 out <<std::endl ;
00361
00362 out << "</DataArray>" <<std::endl ;
00363 out << "</PointData>" <<std::endl ;
00364 out << "</Piece>" <<std::endl ;
00365 out << "</PolyData>" <<std::endl ;
00366 out << "</VTKFile>" <<std::endl ;
00367
00368 QC_TRACE_OUT("END QCPartitioner<TPSystem>::writeInFileVTK");
00369 }
00370
00371 template <class TPSystem>
00372 void QCPartitioner<TPSystem>::writeSubDomainVTK(std::ofstream & out, const QCSubDomain & domain){
00373 QC_TRACE_OUT("BEGIN QCPartitioner<TPSystem>::writeSubDomainVTK");
00374
00375 int nbAtoms = domain.getNbAtoms();
00376
00377 std::vector<int> scalar(nbAtoms,3) ;
00378
00379 QCSubDomainZone * zone = domain.getZone();
00380
00381 for( int j = 0; j < nbAtoms ; ++j){
00382
00383 if(zone[j] == QC_CORE){
00384 scalar[j] = 1 ;
00385 }
00386 else if (this->partType == QC_DIXON_PART ) {
00387
00388 continue;
00389 }
00390 else if(zone[j] == QC_SHELL1 ) {
00391 scalar[j] = 2 ;
00392 }
00393 }
00394 out << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<< std::endl
00395 <<"<PolyData>" << std::endl <<"<Piece NumberOfPoints=\"" << nbAtoms <<"\">"<< std::endl ;
00396 out << "<Points>" <<std::endl ;
00397 out << "<DataArray type=\"Float64\" Name=\"coord_atoms\" NumberOfComponents=\"3\" format=\"ascii\">" <<std::endl ;
00398 QCPoint3D xyz ;
00399 for ( int i = 0 ; i < nbAtoms ; ++i){
00400 domain.getPointAt(i,xyz);
00401 out << xyz[0] << " " <<xyz[1] << " " <<xyz[2] << " " <<std::endl;
00402 }
00403 out << "</DataArray>" <<std::endl ;
00404 out << "</Points>" <<std::endl ;
00405
00406 out << "<PointData Scalars=\"PartitionNumber\"> "<<std::endl ;
00407 out << "<DataArray type=\"Int32\" Name=\"PartitionNumber\" format=\"ascii\"> " <<std::endl;
00408
00409 for ( unsigned int i = 0 ; i < scalar.size() ; ++i){
00410 out << scalar[i] << " " ;
00411 }
00412 out <<std::endl ;
00413
00414 out << "</DataArray>" <<std::endl ;
00415 out << "</PointData>" <<std::endl ;
00416 out << "</Piece>" <<std::endl ;
00417 out << "</PolyData>" <<std::endl ;
00418 out << "</VTKFile>" <<std::endl ;
00419
00420 QC_TRACE_OUT("END QCPartitioner<TPSystem>::writeSubDomainVTK");
00421 }
00422
00426 template class QCPartitioner<TQCMDSystem>;
00427
00428
00429 #if defined(HAVE_MPI) && defined(WITH_MPI_SUPPORT)
00430
00431 #include "QCDistMDSystem.hpp"
00432 template class QCPartitioner<TQCDistMDSystem>;
00433
00434 #endif