ESyS-Particle
4.0.1
|
00001 00002 // // 00003 // Copyright (c) 2003-2011 by The University of Queensland // 00004 // Earth Systems Science Computational Centre (ESSCC) // 00005 // http://www.uq.edu.au/esscc // 00006 // // 00007 // Primary Business: Brisbane, Queensland, Australia // 00008 // Licensed under the Open Software License version 3.0 // 00009 // http://www.opensource.org/licenses/osl-3.0.php // 00010 // // 00012 00013 // -- project includes - 00014 00015 #include "Parallel/SubLattice.h" 00016 #include "Parallel/MpiInfo.h" 00017 #include "Parallel/mpivbuf.h" 00018 #include "Parallel/mpisgvbuf.h" 00019 #include "Parallel/mpibarrier.h" 00020 #include "Parallel/mpia2abuf.h" 00021 #include "Model/BondedInteraction.h" 00022 #include "Model/CappedBondedInteraction.h" 00023 #include "Model/ShortBondedInteraction.h" 00024 #include "Model/DampingIGP.h" 00025 #include "Model/Damping.h" 00026 #include "Model/RotDamping.h" 00027 #include "Model/ABCDampingIGP.h" 00028 #include "Model/ABCDamping.h" 00029 #include "Model/LocalDampingIGP.h" 00030 #include "Model/LocalDamping.h" 00031 #include "Model/RotLocalDamping.h" 00032 #include "Model/FrictionInteraction.h" 00033 #include "Model/FractalFriction.h" 00034 #include "Model/AdhesiveFriction.h" 00035 #include "Model/VWFrictionInteraction.h" 00036 #include "Model/RotFricInteraction.h" 00037 #include "Model/RotElasticInteraction.h" 00038 #include "Model/RotThermFricInteraction.h" 00039 #include "Model/RotThermElasticInteraction.h" 00040 #include "Model/RotThermBondedInteraction.h" 00041 #include "Model/HertzianElasticInteraction.h" 00042 #include "Model/HertzianViscoElasticFrictionInteraction.h" 00043 #include "Model/HertzianViscoElasticInteraction.h" 00044 #include "Model/LinearDashpotInteraction.h" 00045 #include "Model/MeshData.h" 00046 #include "Model/MeshData2D.h" 00047 #include "Model/ETriMeshInteraction.h" 00048 #include "Model/BTriMeshInteraction.h" 00049 #include "Model/BMesh2DInteraction.h" 00050 #include "Model/EMesh2DInteraction.h" 00051 00052 // --- parallel storage includes --- 00053 #include "ppa/src/pp_array.h" 00054 #include "pis/pi_storage_eb.h" 00055 #include "pis/pi_storage_ed.h" 00056 #include "pis/pi_storage_ed_t.h" 00057 #include "pis/pi_storage_ne.h" 00058 #include "pis/pi_storage_ne_t.h" 00059 #include "pis/pi_storage_single.h" 00060 #include "pis/trimesh_pis.h" 00061 #include "pis/trimesh_pis_ne.h" 00062 #include "pis/trimesh_pis_eb.h" 00063 #include "pis/mesh2d_pis_eb.h" 00064 #include "pis/mesh2d_pis_ne.h" 00065 00066 // --- field includes --- 00067 #include "Fields/ScalarParticleFieldSlaveTagged.h" 00068 #include "Fields/VectorParticleFieldSlaveTagged.h" 00069 #include "Fields/ScalarInteractionFieldSlaveTagged.h" 00070 #include "Fields/ScalarParticleFieldSlaveTagged.h" 00071 #include "Fields/ScalarInteractionFieldSlaveTagged.h" 00072 #include "Fields/CheckedScalarInteractionFieldSlave.h" 00073 #include "Fields/CheckedScalarInteractionFieldSlaveTagged.h" 00074 #include "Fields/VectorTriangleFieldSlave.h" 00075 #include "Fields/ScalarTriangleFieldSlave.h" 00076 #include "Fields/VectorEdge2DFieldSlave.h" 00077 00078 #include "Model/BodyForceGroup.h" 00079 00080 #include <mpi.h> 00081 #include <stdlib.h> 00082 #include <stdexcept> 00083 00084 // -- STL includes -- 00085 #include <algorithm> 00086 #include <stdexcept> 00087 #include <boost/limits.hpp> 00088 using std::runtime_error; 00089 00090 // -- IO includes -- 00091 #include <iostream> 00092 using std::cerr; 00093 using std::flush; 00094 using std::endl; 00095 using esys::lsm::CLatticeParam; 00096 00097 //---------------------------------------------- 00098 // TSubLattice functions 00099 //---------------------------------------------- 00100 00108 template <class T> 00109 TSubLattice<T>::TSubLattice( 00110 const CLatticeParam ¶m, 00111 int rank, 00112 MPI_Comm comm, 00113 MPI_Comm worker_comm 00114 ) 00115 : m_ppa(NULL), 00116 m_dpis(), 00117 m_bpis(), 00118 m_singleParticleInteractions(), 00119 m_damping(), 00120 m_WIG(), 00121 m_mesh(), 00122 m_mesh2d(), 00123 m_dt(0), 00124 m_nrange(0), 00125 m_alpha(0), 00126 m_last_ns(0), 00127 m_temp_conn(), 00128 m_rank(0), 00129 m_comm(MPI_COMM_NULL), 00130 m_tml_comm(MPI_COMM_NULL), 00131 m_worker_comm(MPI_COMM_NULL), 00132 m_tml_worker_comm(MPI_COMM_NULL), 00133 m_dims(3, 0), 00134 packtime(0), 00135 unpacktime(0), 00136 commtime(0.0), 00137 forcetime(0.0), 00138 m_field_slaves(), 00139 m_pTimers(NULL) 00140 { 00141 // cout << "TSubLattice<T>::TSubLattice at " << rank << endl << flush; 00142 // -- MPI stuff -- 00143 m_rank=rank; 00144 00145 // set global communicator 00146 m_comm=comm; 00147 m_tml_comm.setComm(m_comm); 00148 00149 m_dims = param.processDims(); 00150 00151 m_worker_comm=worker_comm; 00152 // MPI_Comm_size(m_worker_comm,&m_num_workers); 00153 m_tml_worker_comm.setComm(m_worker_comm); 00154 00155 00156 // -- set parameters 00157 m_alpha=param.alpha(); 00158 m_nrange=param.nrange(); 00159 // cout << "dt,nrange,alpha : " << m_dt << " , " << m_nrange << " , " << m_alpha << "\n"; 00160 00161 commtime=0.0; 00162 packtime=0.0; 00163 unpacktime=0.0; 00164 forcetime=0.0; 00165 00166 m_last_ns=-1; 00167 } 00168 00169 template <class T> 00170 TSubLattice<T>::~TSubLattice() 00171 { 00172 console.Debug() << "TSubLattice<T>::~TSubLattice(): enter\n"; 00173 console.Debug() 00174 << "TSubLattice<T>::~TSubLattice():" 00175 << " deleting wall interaction groups...\n"; 00176 for( 00177 typename map<string,AWallInteractionGroup<T>*>::iterator vit=m_WIG.begin(); 00178 vit!=m_WIG.end(); 00179 vit++ 00180 ) 00181 { 00182 delete vit->second; 00183 } 00184 console.Debug() 00185 << "TSubLattice<T>::~TSubLattice():" 00186 << " deleting particle array...\n"; 00187 if (m_ppa != NULL) delete m_ppa; 00188 console.Debug() << "TSubLattice<T>::~TSubLattice(): exit\n"; 00189 } 00190 00197 template <class T> 00198 void TSubLattice<T>::initNeighborTable(const Vec3& min,const Vec3& max) 00199 { 00200 console.XDebug() << "TSubLattice<T>::initNeighborTable(" << min << "," << max << ")\n"; 00201 // make size fit range 00202 double xsize=max.X()-min.X(); 00203 xsize=m_nrange*ceil(xsize/m_nrange); 00204 double ysize=max.Y()-min.Y(); 00205 ysize=m_nrange*ceil(ysize/m_nrange); 00206 double zsize=max.Z()-min.Z(); 00207 zsize=m_nrange*ceil(zsize/m_nrange); 00208 Vec3 grow=Vec3(xsize,ysize,zsize)-(max-min); // size increase 00209 Vec3 nmin=min-0.5*grow; // distribute symmetrically 00210 Vec3 nmax=max+0.5*grow; 00211 console.XDebug() << "range=" << m_nrange << ", new min,max: " << nmin << ", " << nmax << "\n"; 00212 00213 // construct particle array 00214 TML_Comm *ntcomm=new TML_Comm(m_worker_comm); 00215 m_ppa=new ParallelParticleArray<T>(ntcomm,m_dims,nmin,nmax,m_nrange,m_alpha); 00216 //m_ppa=new ParallelParticleArray<T>(ntcomm,3,nmin,nmax,m_nrange); 00217 00218 // makeFields(); // put here to make sure ppa is initialized before makeFields 00219 00220 console.XDebug() << "end TSubLattice<T>::initNeighborTable\n"; 00221 } 00222 00223 template <class T> 00224 void TSubLattice<T>::do2dCalculations(bool do2d) 00225 { 00226 T::setDo2dCalculations(do2d); 00227 } 00228 00229 template <class T> 00230 int TSubLattice<T>::getNumParticles() 00231 { 00232 return m_ppa->getInnerSize(); 00233 } 00234 00242 template <class T> 00243 void TSubLattice<T>::initNeighborTable(const Vec3& min,const Vec3& max,const vector<bool>& circ) 00244 { 00245 console.XDebug() << "TSubLattice<T>::initNeighborTable(" << min << "," << max << ") circ\n"; 00246 double xsize,ysize,zsize; 00247 // if dimension is circular, check if size fits range, otherwise make it fit 00248 // x - dim 00249 if(circ[0]) 00250 { 00251 xsize=max.X()-min.X(); 00252 if(fabs(xsize/m_nrange-lrint(xsize/m_nrange))>1e-6) 00253 { 00254 //console.Critical() << "circular x-dimension doesn't fit range !\n"; 00255 console.Info() << "Circular x-size incompatible with range, adjusting...\n"; 00256 m_nrange = xsize/floor(xsize/m_nrange); 00257 console.Info() << "New range = " << m_nrange << "\n"; 00258 } 00259 //xsize+=2.0*m_nrange; // padding on the circular ends 00260 } 00261 else 00262 { 00263 xsize=max.X()-min.X(); 00264 xsize=m_nrange*ceil(xsize/m_nrange); 00265 } 00266 // y - dim 00267 if(circ[1]) 00268 { 00269 ysize=max.Y()-min.Y(); 00270 if(fabs(ysize/m_nrange-lrint(ysize/m_nrange))>1e-6) 00271 { 00272 console.Critical() << "circular y-dimension doesn't fit range !\n"; 00273 } 00274 ysize+=2.0*m_nrange; // padding on the circular ends 00275 } 00276 else 00277 { 00278 ysize=max.Y()-min.Y(); 00279 ysize=m_nrange*ceil(ysize/m_nrange); 00280 } 00281 // z - dim 00282 if(circ[2]) 00283 { 00284 zsize=max.Z()-min.Z(); 00285 if(fabs(zsize/m_nrange-lrint(zsize/m_nrange))>1e-6) 00286 { 00287 console.Critical() << "circular z-dimension doesn't fit range !\n"; 00288 } 00289 zsize+=2.0*m_nrange; // padding on the circular ends 00290 } 00291 else 00292 { 00293 zsize=max.Z()-min.Z(); 00294 zsize=m_nrange*ceil(zsize/m_nrange); 00295 } 00296 Vec3 grow=Vec3(xsize,ysize,zsize)-(max-min); // size increase 00297 Vec3 nmin=min-0.5*grow; // distribute symmetrically 00298 Vec3 nmax=max+0.5*grow; 00299 console.XDebug() << "range, new min, max: " << m_nrange << " " << nmin << nmax << "\n"; 00300 // construct particle array 00301 TML_Comm *ntcomm=new TML_Comm(m_worker_comm); 00302 m_ppa=new ParallelParticleArray<T>(ntcomm,m_dims,circ,nmin,nmax,m_nrange,m_alpha); 00303 00304 // makeFields(); // put here to make sure ppa is initialized before makeFields 00305 00306 console.XDebug() << "end TSubLattice<T>::initNeighborTable (circ)\n"; 00307 } 00308 00315 template <class T> 00316 void TSubLattice<T>::receiveParticles() 00317 { 00318 console.XDebug() << "TSubLattice<T>::receiveParticles: enter\n"; 00319 00320 vector<T> recv_buffer; 00321 CMPIBarrier barrier(m_comm); 00322 00323 m_tml_comm.recv_broadcast_cont_packed(recv_buffer,0); 00324 console.XDebug() << "recvd " << recv_buffer.size() << " particles \n"; 00325 m_ppa->insert(recv_buffer); 00326 00327 barrier.wait("TSubLattice<T>::receiveParticles"); 00328 00329 console.XDebug() << "TSubLattice<T>::receiveParticles: exit\n"; 00330 } 00331 00332 00338 template <class T> 00339 void TSubLattice<T>::receiveConnections() 00340 { 00341 console.XDebug() << "TSubLattice<T>::receiveConnections: enter\n"; 00342 00343 vector<int> recv_buffer; 00344 CMPIBarrier barrier(m_comm); 00345 00346 m_tml_comm.recv_broadcast_cont_packed(recv_buffer,0); 00347 console.XDebug() << "recvd " << recv_buffer.size() << " connections \n"; 00348 vector<int>::iterator it; 00349 for (it = recv_buffer.begin(); it != recv_buffer.end(); it+=3) 00350 { 00351 if ( (m_ppa->getParticlePtrByIndex( *(it+1)) == NULL ) || 00352 (m_ppa->getParticlePtrByIndex( *(it+2)) == NULL ) ) 00353 { 00354 continue; 00355 } 00356 m_temp_conn[*(it)].push_back(*(it+1)); 00357 m_temp_conn[*(it)].push_back(*(it+2)); 00358 } 00359 00360 barrier.wait("TSubLattice<T>::receiveConnections"); 00361 00362 console.XDebug() << "TSubLattice<T>::receiveConnections: exit\n"; 00363 } 00364 00365 00369 template <class T> 00370 void TSubLattice<T>::addWall() 00371 { 00372 console.XDebug() << "TSubLattice<T>::addWall: enter\n" ; 00373 CVarMPIBuffer param_buffer(m_comm); 00374 param_buffer.receiveBroadcast(0); 00375 00376 string name=param_buffer.pop_string(); 00377 Vec3 ipos=param_buffer.pop_vector(); 00378 Vec3 inorm=param_buffer.pop_vector(); 00379 00380 m_walls[name]=new CWall(ipos,inorm); 00381 00382 console.XDebug() << "TSubLattice<T>::addWall: exit\n" ; 00383 } 00384 00388 template <class T> 00389 void TSubLattice<T>::addElasticWIG() 00390 { 00391 console.XDebug() << "TSubLattice<T>::addElasticWIG: enter\n" ; 00392 CVarMPIBuffer param_buffer(m_comm); 00393 00394 // receive params from master 00395 param_buffer.receiveBroadcast(0); 00396 00397 CEWallIGP* wigp=extractEWallIGP(¶m_buffer); 00398 00399 string wallname=wigp->getWallName(); 00400 map<string,CWall*>::iterator iter=m_walls.find(wallname); 00401 if(iter!=m_walls.end()){ 00402 AWallInteractionGroup<T>* newCEWIG = 00403 new CEWallInteractionGroup<T>( 00404 &m_tml_worker_comm, 00405 m_walls[wallname], 00406 wigp 00407 ); 00408 newCEWIG->Update(m_ppa); 00409 m_WIG.insert(make_pair(wigp->getName(),newCEWIG)); 00410 } else { 00411 std::stringstream msg; 00412 msg << "wall name '" << wallname << "' not found in map of walls"; 00413 throw std::runtime_error(msg.str().c_str()); 00414 } 00415 00416 delete wigp; 00417 console.XDebug() << "TSubLattice<T>::addElasticWIG: exit\n" ; 00418 } 00419 00420 00424 template <class T> 00425 void TSubLattice<T>::addBondedWIG() 00426 { 00427 console.XDebug() << "TSubLattice<T>::addBondedWIG: enter\n" ; 00428 CVarMPIBuffer param_buffer(m_comm); 00429 00430 // receive params from master 00431 param_buffer.receiveBroadcast(0); 00432 00433 CBWallIGP* wigp=extractBWallIGP(¶m_buffer); 00434 00435 string wallname=wigp->getWallName(); 00436 map<string,CWall*>::iterator iter=m_walls.find(wallname); 00437 if(iter!=m_walls.end()){ 00438 AWallInteractionGroup<T>* newCBWIG=new CBWallInteractionGroup<T>(&m_tml_worker_comm,m_walls[wallname],wigp); 00439 newCBWIG->Update(m_ppa); 00440 m_WIG.insert(make_pair(wigp->getName(),newCBWIG)); 00441 } else { 00442 console.Error() << "wall name " << wallname << " not found in map of walls\n"; 00443 } 00444 00445 delete wigp; 00446 console.XDebug() << "TSubLattice<T>::addBondedWIG: exit\n" ; 00447 } 00448 00452 template <class T> 00453 void TSubLattice<T>::addDirBondedWIG() 00454 { 00455 console.XDebug() << "TSubLattice<T>::addDirBondedWIG: enter\n" ; 00456 CVarMPIBuffer param_buffer(m_comm); 00457 00458 // receive params from master 00459 param_buffer.receiveBroadcast(0); 00460 00461 CSoftBWallIGP* wigp=extractSoftBWallIGP(¶m_buffer); 00462 00463 string wallname=wigp->getWallName(); 00464 map<string,CWall*>::iterator iter=m_walls.find(wallname); 00465 if(iter!=m_walls.end()){ 00466 AWallInteractionGroup<T>* newCDWIG=new CSoftBWallInteractionGroup<T>(&m_tml_worker_comm,m_walls[wallname],wigp); 00467 newCDWIG->Update(m_ppa); 00468 m_WIG.insert(make_pair(wigp->getName(),newCDWIG)); 00469 } else { 00470 console.Error() << "wall name " << wallname << " not found in map of walls\n"; 00471 } 00472 00473 delete wigp; 00474 console.XDebug() << "TSubLattice<T>::addDirBondedWIG: exit\n" ; 00475 } 00476 00480 template <class T> 00481 void TSubLattice<T>::getWallPos() 00482 { 00483 console.XDebug() << "TSubLattice<T>::getWallPosition: enter\n" ; 00484 CVarMPIBuffer param_buffer(m_comm); 00485 Vec3 pos; 00486 00487 // receive params from master 00488 param_buffer.receiveBroadcast(0); 00489 00490 std::string wname=param_buffer.pop_string(); 00491 console.XDebug() << "Wall name: " << wname << "\n"; 00492 00493 // find wall 00494 map<string,CWall*>::iterator iter=m_walls.find(wname); 00495 if(iter!=m_walls.end()){ 00496 pos=(iter->second)->getPos(); 00497 console.XDebug() << "Wall pos: " << pos << "\n"; 00498 } else { 00499 pos=Vec3(0.0,0.0,0.0); 00500 } 00501 00502 vector<Vec3> vpos; 00503 vpos.push_back(pos); 00504 m_tml_comm.send_gather(vpos,0); 00505 console.XDebug() << "TSubLattice<T>::getWallPosition: exit\n" ; 00506 } 00507 00511 template <class T> 00512 void TSubLattice<T>::getWallForce() 00513 { 00514 console.XDebug() << "TSubLattice<T>::getWallForce: enter\n" ; 00515 CVarMPIBuffer param_buffer(m_comm); 00516 Vec3 force; 00517 00518 // receive params from master 00519 param_buffer.receiveBroadcast(0); 00520 00521 std::string wname=param_buffer.pop_string(); 00522 console.XDebug() << "Wall name: " << wname << "\n"; 00523 00524 // find wall 00525 map<string,CWall*>::iterator iter=m_walls.find(wname); 00526 if(iter!=m_walls.end()){ 00527 force=(iter->second)->getForce(); 00528 console.XDebug() << "Wall force: " << force << "\n"; 00529 } else { 00530 force=Vec3(0.0,0.0,0.0); 00531 } 00532 00533 vector<Vec3> vforce; 00534 vforce.push_back(force); 00535 m_tml_comm.send_gather(vforce,0); 00536 console.XDebug() << "TSubLattice<T>::getWallForce: exit\n" ; 00537 } 00538 00542 template <class T> 00543 void TSubLattice<T>::addViscWIG() 00544 { 00545 console.XDebug() << "TSubLattice<T>::addViscWIG: enter\n" ; 00546 CVarMPIBuffer param_buffer(m_comm); 00547 00548 // receive params from master 00549 param_buffer.receiveBroadcast(0); 00550 00551 CVWallIGP* wigp=extractVWallIGP(¶m_buffer); 00552 00553 string wallname=wigp->getWallName(); 00554 map<string,CWall*>::iterator iter=m_walls.find(wallname); 00555 if(iter!=m_walls.end()){ 00556 AWallInteractionGroup<T>* newCVWIG=new CViscWallIG<T>(&m_tml_worker_comm,m_walls[wallname],wigp); 00557 newCVWIG->Update(m_ppa); 00558 m_WIG.insert(make_pair(wigp->getName(),newCVWIG)); 00559 } else { 00560 console.Error() << "wall name " << wallname << " not found in map of walls\n"; 00561 } 00562 00563 delete wigp; 00564 console.XDebug() << "TSubLattice<T>::addViscWIG: exit\n" ; 00565 } 00566 00570 template <class T> 00571 void TSubLattice<T>::addPairIG() 00572 { 00573 console.XDebug() << "TSubLattice<T>::addPairIG()\n"; 00574 CVarMPIBuffer param_buffer(m_comm,2000); 00575 00576 // get params 00577 param_buffer.receiveBroadcast(0); 00578 string type = param_buffer.pop_string(); 00579 console.XDebug()<< "PIG type: " << type.c_str() << "\n"; 00580 string name = param_buffer.pop_string(); 00581 console.XDebug()<< "PIG name: " << name.c_str() << "\n"; 00582 00583 doAddPIG(name,type,param_buffer,false); 00584 00585 console.XDebug() << "end TSubLattice<T>::addPairIG()\n"; 00586 } 00587 00591 template <class T> 00592 void TSubLattice<T>::addTaggedPairIG() 00593 { 00594 console.XDebug() << "TSubLattice<T>::addTaggedPairIG()\n"; 00595 CVarMPIBuffer param_buffer(m_comm,2000); 00596 00597 // get params 00598 param_buffer.receiveBroadcast(0); 00599 string type = param_buffer.pop_string(); 00600 console.XDebug()<< "PIG type: " << type.c_str() << "\n"; 00601 string name = param_buffer.pop_string(); 00602 console.XDebug()<< "PIG name: " << name.c_str() << "\n"; 00603 00604 doAddPIG(name,type,param_buffer,true); 00605 00606 console.XDebug() << "end TSubLattice<T>::addTaggedPairIG()\n"; 00607 } 00608 00616 template <class T> 00617 bool TSubLattice<T>::doAddPIG(const string& name,const string& type,CVarMPIBuffer& param_buffer, bool tagged) 00618 { 00619 bool res=false; 00620 AParallelInteractionStorage* new_pis = NULL; 00621 00622 if(type=="Elastic") { 00623 CElasticIGP eigp; 00624 eigp.m_k=param_buffer.pop_double(); 00625 eigp.m_scaling=static_cast<bool>(param_buffer.pop_int()); 00626 if(tagged){ 00627 int tag1=param_buffer.pop_int(); 00628 int mask1=param_buffer.pop_int(); 00629 int tag2=param_buffer.pop_int(); 00630 int mask2=param_buffer.pop_int(); 00631 console.XDebug() << "tag1, mask1, tag2, mask2 " 00632 << tag1 << " , " << mask1 << " , " 00633 << tag2 << " , " << mask2 << "\n"; 00634 new_pis=new ParallelInteractionStorage_NE_T<T,CElasticInteraction>(m_ppa,eigp,tag1,mask1,tag2,mask2); 00635 }else{ 00636 new_pis=new ParallelInteractionStorage_NE<T,CElasticInteraction>(m_ppa,eigp); 00637 } 00638 res=true; 00639 } else if (type=="Friction") { 00640 CFrictionIGP figp; 00641 figp.k=param_buffer.pop_double(); 00642 figp.mu=param_buffer.pop_double(); 00643 figp.k_s=param_buffer.pop_double(); 00644 figp.dt=param_buffer.pop_double(); 00645 figp.m_scaling=static_cast<bool>(param_buffer.pop_int()); 00646 console.XDebug() << "k,mu,k_s,dt: " << figp.k << " , " << figp.mu << " , " 00647 << figp.k_s << " , " << figp.dt << "\n"; 00648 new_pis=new ParallelInteractionStorage_ED<T,CFrictionInteraction>(m_ppa,figp); 00649 res=true; 00650 } else if (type=="AdhesiveFriction") { 00651 CAdhesiveFrictionIGP figp; 00652 figp.k=param_buffer.pop_double(); 00653 figp.mu=param_buffer.pop_double(); 00654 figp.k_s=param_buffer.pop_double(); 00655 figp.dt=param_buffer.pop_double(); 00656 figp.r_cut=param_buffer.pop_double(); 00657 console.XDebug() 00658 << "k,mu,k_s,dt,r_cut: " << figp.k << " , " << figp.mu << " , " 00659 << figp.k_s << " , " << figp.dt << " " << figp.r_cut << "\n"; 00660 new_pis=new ParallelInteractionStorage_ED<T,CAdhesiveFriction>(m_ppa,figp); 00661 res=true; 00662 } else if (type=="FractalFriction") { 00663 FractalFrictionIGP figp; 00664 figp.k=param_buffer.pop_double(); 00665 figp.mu_0=param_buffer.pop_double(); 00666 figp.k_s=param_buffer.pop_double(); 00667 figp.dt=param_buffer.pop_double(); 00668 console.XDebug() << "k,mu_0,k_s,dt: " << figp.k << " , " << figp.mu_0 << " , " 00669 << figp.k_s << " , " << figp.dt << "\n"; 00670 figp.x0=param_buffer.pop_double(); 00671 figp.y0=param_buffer.pop_double(); 00672 figp.dx=param_buffer.pop_double(); 00673 figp.dy=param_buffer.pop_double(); 00674 figp.nx=param_buffer.pop_int(); 00675 figp.ny=param_buffer.pop_int(); 00676 console.XDebug() 00677 <<"x0,y0,dx,dy,nx,ny: " 00678 << figp.x0 << " , " << figp.y0 << " , " 00679 << figp.dx << " , " << figp.dy << " ," 00680 << figp.nx << " , " << figp.ny << "\n"; 00681 figp.mu = boost::shared_ptr<double>(new double[figp.nx*figp.ny]); 00682 00683 for(int i=0;i<figp.nx*figp.ny;i++) 00684 { 00685 (figp.mu.get())[i]=param_buffer.pop_double(); 00686 // console.XDebug() << i << " , " << figp.mu[i] << "\n"; 00687 } 00688 new_pis = new ParallelInteractionStorage_ED<T,CFractalFriction>(m_ppa,figp); 00689 res=true; 00690 } else if(type=="VWFriction") { 00691 VWFrictionIGP figp; 00692 00693 figp.k=param_buffer.pop_double(); 00694 figp.mu=param_buffer.pop_double(); 00695 figp.k_s=param_buffer.pop_double(); 00696 figp.dt=param_buffer.pop_double(); 00697 figp.m_alpha=param_buffer.pop_double(); 00698 console.XDebug() 00699 << "k,mu,k_s,dt,alpha: " << figp.k << " , " << figp.mu << " , " 00700 << figp.k_s << " , " << figp.dt << "\n"; 00701 new_pis=new ParallelInteractionStorage_ED<T,CVWFriction>(m_ppa,figp); 00702 res=true; 00703 } else if(type=="RotElastic"){ 00704 CRotElasticIGP reigp; 00705 reigp.m_kr=param_buffer.pop_double(); 00706 new_pis=new ParallelInteractionStorage_NE<CRotParticle,CRotElasticInteraction>(m_ppa,reigp); 00707 res=true; 00708 } else if (type=="RotFriction"){ 00709 CRotFrictionIGP rfigp; 00710 rfigp.k=param_buffer.pop_double(); 00711 rfigp.mu_s=param_buffer.pop_double(); 00712 rfigp.mu_d=param_buffer.pop_double(); 00713 rfigp.k_s=param_buffer.pop_double(); 00714 rfigp.dt=param_buffer.pop_double(); 00715 rfigp.scaling=static_cast<bool>(param_buffer.pop_int()); 00716 console.XDebug() 00717 << "k,mu_s,mu_d,k_s,dt,scaling: " << rfigp.k << " , " 00718 << rfigp.mu_s << " , " << rfigp.mu_d << " , " 00719 << rfigp.k_s << " , " << rfigp.dt << " , " << rfigp.scaling << "\n"; 00720 if(tagged){ 00721 int tag1=param_buffer.pop_int(); 00722 int mask1=param_buffer.pop_int(); 00723 int tag2=param_buffer.pop_int(); 00724 int mask2=param_buffer.pop_int(); 00725 console.XDebug() << "tag1, mask1, tag2, mask2 " 00726 << tag1 << " , " << mask1 << " , " 00727 << tag2 << " , " << mask2 << "\n"; 00728 new_pis=new ParallelInteractionStorage_ED_T<CRotParticle,CRotFrictionInteraction>(m_ppa,rfigp,tag1,mask1,tag2,mask2); 00729 } else { 00730 new_pis=new ParallelInteractionStorage_ED<CRotParticle,CRotFrictionInteraction>(m_ppa,rfigp); 00731 } 00732 res=true; 00733 } else if (type == "RotThermElastic") { 00734 CRotThermElasticIGP eigp; 00735 eigp.m_kr = param_buffer.pop_double(); 00736 eigp.diffusivity = param_buffer.pop_double(); 00737 console.XDebug() 00738 << "k=" << eigp.m_kr << " , " 00739 << "diffusivity=" << eigp.diffusivity << "\n"; 00740 00741 new_pis = 00742 new ParallelInteractionStorage_NE< 00743 CRotThermParticle, 00744 CRotThermElasticInteraction 00745 >( 00746 m_ppa, 00747 eigp 00748 ); 00749 res=true; 00750 } else if (type == "RotThermFriction") { 00751 CRotThermFrictionIGP rfigp; 00752 rfigp.k=param_buffer.pop_double(); 00753 rfigp.mu_s=param_buffer.pop_double(); 00754 rfigp.mu_d=param_buffer.pop_double(); 00755 rfigp.k_s=param_buffer.pop_double(); 00756 rfigp.diffusivity=param_buffer.pop_double(); 00757 rfigp.dt=param_buffer.pop_double(); 00758 console.XDebug() 00759 << "k=" << rfigp.k << " , " 00760 << "mu_d=" << rfigp.mu_d << " , " 00761 << "mu_s=" << rfigp.mu_s << " , " 00762 << "k_s=" << rfigp.k_s << " , " 00763 << "diffusivity=" << rfigp.diffusivity << " , " 00764 << "dt=" << rfigp.dt << "\n"; 00765 00766 new_pis = 00767 new ParallelInteractionStorage_ED< 00768 CRotThermParticle, 00769 CRotThermFrictionInteraction 00770 >( 00771 m_ppa, 00772 rfigp 00773 ); 00774 res=true; 00775 } else if(type=="HertzianElastic") { 00776 CHertzianElasticIGP heigp; 00777 heigp.m_E=param_buffer.pop_double(); 00778 heigp.m_nu=param_buffer.pop_double(); 00779 if(tagged){ 00780 int tag1=param_buffer.pop_int(); 00781 int mask1=param_buffer.pop_int(); 00782 int tag2=param_buffer.pop_int(); 00783 int mask2=param_buffer.pop_int(); 00784 new_pis=new ParallelInteractionStorage_NE_T<T,CHertzianElasticInteraction>(m_ppa,heigp,tag1,mask1,tag2,mask2); 00785 }else{ 00786 new_pis=new ParallelInteractionStorage_NE<T,CHertzianElasticInteraction>(m_ppa,heigp); 00787 } 00788 res=true; 00789 } else if(type=="HertzianViscoElasticFriction") { 00790 CHertzianViscoElasticFrictionIGP hvefigp; 00791 hvefigp.m_A=param_buffer.pop_double(); 00792 hvefigp.m_E=param_buffer.pop_double(); 00793 hvefigp.m_nu=param_buffer.pop_double(); 00794 hvefigp.mu=param_buffer.pop_double(); 00795 hvefigp.k_s=param_buffer.pop_double(); 00796 hvefigp.dt=param_buffer.pop_double(); 00797 if(tagged){ 00798 int tag1=param_buffer.pop_int(); 00799 int mask1=param_buffer.pop_int(); 00800 int tag2=param_buffer.pop_int(); 00801 int mask2=param_buffer.pop_int(); 00802 new_pis=new ParallelInteractionStorage_NE_T<T,CHertzianViscoElasticFrictionInteraction>(m_ppa,hvefigp,tag1,mask1,tag2,mask2); 00803 }else{ 00804 new_pis=new ParallelInteractionStorage_NE<T,CHertzianViscoElasticFrictionInteraction>(m_ppa,hvefigp); 00805 } 00806 res=true; 00807 } else if(type=="HertzianViscoElastic") { 00808 CHertzianViscoElasticIGP hveigp; 00809 hveigp.m_A=param_buffer.pop_double(); 00810 hveigp.m_E=param_buffer.pop_double(); 00811 hveigp.m_nu=param_buffer.pop_double(); 00812 if(tagged){ 00813 int tag1=param_buffer.pop_int(); 00814 int mask1=param_buffer.pop_int(); 00815 int tag2=param_buffer.pop_int(); 00816 int mask2=param_buffer.pop_int(); 00817 new_pis=new ParallelInteractionStorage_NE_T<T,CHertzianViscoElasticInteraction>(m_ppa,hveigp,tag1,mask1,tag2,mask2); 00818 }else{ 00819 new_pis=new ParallelInteractionStorage_NE<T,CHertzianViscoElasticInteraction>(m_ppa,hveigp); 00820 } 00821 res=true; 00822 } else if(type=="LinearDashpot") { 00823 CLinearDashpotIGP heigp; 00824 heigp.m_damp=param_buffer.pop_double(); 00825 heigp.m_cutoff=param_buffer.pop_double(); 00826 if(tagged){ 00827 int tag1=param_buffer.pop_int(); 00828 int mask1=param_buffer.pop_int(); 00829 int tag2=param_buffer.pop_int(); 00830 int mask2=param_buffer.pop_int(); 00831 new_pis=new ParallelInteractionStorage_NE_T<T,CLinearDashpotInteraction>(m_ppa,heigp,tag1,mask1,tag2,mask2); 00832 }else{ 00833 new_pis=new ParallelInteractionStorage_NE<T,CLinearDashpotInteraction>(m_ppa,heigp); 00834 } 00835 res=true; 00836 } else { 00837 cerr << "Unknown interaction group name " 00838 << type 00839 << " in TSubLattice<T>::addPairIG()" << endl; 00840 } 00841 00842 // add InteractionGroup to map 00843 if(res) m_dpis.insert(make_pair(name,new_pis)); 00844 00845 return res; 00846 } 00847 00848 00852 template <class T> 00853 void TSubLattice<T>::addTriMesh() 00854 { 00855 console.XDebug() << "TSubLattice<T>::addTriMesh()\n"; 00856 00857 // MPI buffer 00858 CVarMPIBuffer param_buffer(m_comm); 00859 // data buffers 00860 vector<MeshNodeData> node_recv_buffer; 00861 vector<MeshTriData> tri_recv_buffer; 00862 00863 // receive params 00864 param_buffer.receiveBroadcast(0); 00865 00866 // extract name 00867 string name = param_buffer.pop_string(); 00868 console.XDebug()<< "TriMesh name: " << name.c_str() << "\n"; 00869 00870 // receive nodes 00871 m_tml_comm.recv_broadcast_cont_packed(node_recv_buffer,0); 00872 console.XDebug() << "recvd " << node_recv_buffer.size() << " nodes \n"; 00873 00874 // receive triangles 00875 m_tml_comm.recv_broadcast_cont_packed(tri_recv_buffer,0); 00876 console.XDebug() << "recvd " << tri_recv_buffer.size() << " triangles \n"; 00877 00878 // load mesh into new TriMesh 00879 TriMesh* new_tm=new TriMesh(); 00880 new_tm->LoadMesh(node_recv_buffer,tri_recv_buffer); 00881 00882 m_mesh.insert(make_pair(name,new_tm)); 00883 00884 console.XDebug() << "end TSubLattice<T>::addTriMesh()\n"; 00885 } 00886 00890 template <class T> 00891 void TSubLattice<T>::addTriMeshIG() 00892 { 00893 console.XDebug() << "TSubLattice<T>::addTriMeshIG()\n"; 00894 00895 // MPI buffer 00896 CVarMPIBuffer param_buffer(m_comm); 00897 00898 // receive params 00899 param_buffer.receiveBroadcast(0); 00900 00901 // extract name & type 00902 string type = param_buffer.pop_string(); 00903 console.XDebug()<< "TriMeshIG type: " << type.c_str() << "\n"; 00904 string name = param_buffer.pop_string(); 00905 console.XDebug()<< "TriMeshIG name: " << name.c_str() << "\n"; 00906 string meshname = param_buffer.pop_string(); 00907 console.XDebug()<< "TriMeshIG mesh name: " << meshname.c_str() << "\n"; 00908 00909 // get pointer to mesh 00910 TriMesh* tmp=NULL; 00911 if (m_mesh.find(meshname) != m_mesh.end()) 00912 { 00913 tmp = m_mesh[meshname]; 00914 } 00915 if(tmp==NULL){ 00916 throw runtime_error("unknown mesh name in TSubLattice<T>::addTriMeshIG:" + meshname); 00917 } 00918 // switch on type,extract params & construc new TriMeshIG 00919 if(type=="Elastic") 00920 { 00921 AParallelInteractionStorage* new_pis; 00922 ETriMeshIP tmi; 00923 tmi.k=param_buffer.pop_double(); 00924 new_pis = new TriMesh_PIS_NE<T,ETriMeshInteraction>(tmp,m_ppa,tmi); 00925 m_dpis.insert(make_pair(name,new_pis)); 00926 } else { // unknown type-> throw 00927 throw runtime_error("unknown type in TSubLattice<T>::addTriMeshIG:" + type); 00928 } 00929 00930 00931 console.XDebug() << "end TSubLattice<T>::addTriMeshIG()\n"; 00932 } 00933 00937 template <class T> 00938 void TSubLattice<T>::addBondedTriMeshIG() 00939 { 00940 console.XDebug() << "TSubLattice<T>::addBondedTriMeshIG()\n"; 00941 00942 // MPI buffer 00943 CVarMPIBuffer param_buffer(m_comm); 00944 00945 // receive params 00946 BTriMeshIP param; 00947 param_buffer.receiveBroadcast(0); 00948 00949 // extract name.meshname & params 00950 string name = param_buffer.pop_string(); 00951 console.XDebug()<< "BTriMeshIG name: " << name.c_str() << "\n"; 00952 string meshname = param_buffer.pop_string(); 00953 console.XDebug()<< "BTriMeshIG mesh name: " << meshname.c_str() << "\n"; 00954 param.k = param_buffer.pop_double(); 00955 console.XDebug()<< "BTriMeshIG k : " << param.k << "\n"; 00956 param.brk = param_buffer.pop_double(); 00957 console.XDebug()<< "BTriMeshIG r_break: " << param.brk << "\n"; 00958 string buildtype = param_buffer.pop_string(); 00959 console.XDebug()<< "BTriMeshIG build type: " << buildtype.c_str() << "\n"; 00960 00961 // get pointer to mesh 00962 TriMesh* tmp=NULL; 00963 if (m_mesh.find(meshname) != m_mesh.end()) 00964 { 00965 tmp = m_mesh[meshname]; 00966 } 00967 if(tmp==NULL){ 00968 throw runtime_error("unknown mesh name in TSubLattice<T>::addTriMeshIG:" + meshname); 00969 } 00970 00971 // setup new interaction storage 00972 TriMesh_PIS_EB<T,BTriMeshInteraction>* new_pis=new TriMesh_PIS_EB<T,BTriMeshInteraction>(tmp,m_ppa,param); 00973 // switch on buildtype, extract buildparam & build 00974 if(buildtype=="BuildByTag"){ 00975 int tag=param_buffer.pop_int(); 00976 int mask=param_buffer.pop_int(); 00977 new_pis->buildFromPPATagged(tag,mask); 00978 m_bpis.insert(make_pair(name,new_pis)); 00979 } else if(buildtype=="BuildByGap"){ 00980 double max_gap=param_buffer.pop_double(); 00981 new_pis->buildFromPPAByGap(max_gap); 00982 m_bpis.insert(make_pair(name,new_pis)); 00983 } else { // unknown type-> throw 00984 throw runtime_error("unknown build type in TSubLattice<T>::addBondedTriMeshIG:" + buildtype); 00985 } 00986 00987 console.XDebug() << "end TSubLattice<T>::addBondedTriMeshIG()\n"; 00988 } 00989 00993 template <class T> 00994 void TSubLattice<T>::addMesh2D() 00995 { 00996 console.XDebug() << "TSubLattice<T>::addMesh2D()\n"; 00997 00998 // MPI buffer 00999 CVarMPIBuffer param_buffer(m_comm); 01000 // data buffers 01001 vector<MeshNodeData2D> node_recv_buffer; 01002 vector<MeshEdgeData2D> edge_recv_buffer; 01003 01004 // receive params 01005 param_buffer.receiveBroadcast(0); 01006 01007 // extract name 01008 string name = param_buffer.pop_string(); 01009 console.XDebug()<< "Mesh2D name: " << name.c_str() << "\n"; 01010 01011 // receive nodes 01012 m_tml_comm.recv_broadcast_cont_packed(node_recv_buffer,0); 01013 console.XDebug() << "recvd " << node_recv_buffer.size() << " nodes \n"; 01014 01015 // receive edges 01016 m_tml_comm.recv_broadcast_cont_packed(edge_recv_buffer,0); 01017 console.XDebug() << "recvd " << edge_recv_buffer.size() << " edges \n"; 01018 01019 // load mesh into new 2D Mesh 01020 Mesh2D* new_tm=new Mesh2D(); 01021 new_tm->LoadMesh(node_recv_buffer,edge_recv_buffer); 01022 01023 m_mesh2d.insert(make_pair(name,new_tm)); 01024 01025 console.XDebug() << "end TSubLattice<T>::addMesh2D()\n"; 01026 } 01027 01032 template <class T> 01033 void TSubLattice<T>::addMesh2DIG() 01034 { 01035 console.XDebug() << "TSubLattice<T>::addMesh2DIG()\n"; 01036 01037 // MPI buffer 01038 CVarMPIBuffer param_buffer(m_comm); 01039 01040 // receive params 01041 param_buffer.receiveBroadcast(0); 01042 01043 // extract name & type 01044 string type = param_buffer.pop_string(); 01045 console.XDebug()<< "Mesh2DIG type: " << type.c_str() << "\n"; 01046 string name = param_buffer.pop_string(); 01047 console.XDebug()<< "Mesh2DIG name: " << name.c_str() << "\n"; 01048 string meshname = param_buffer.pop_string(); 01049 console.XDebug()<< "Mesh2DIG mesh name: " << meshname.c_str() << "\n"; 01050 01051 // get pointer to mesh 01052 Mesh2D* tmp=NULL; 01053 if (m_mesh2d.find(meshname) != m_mesh2d.end()) 01054 { 01055 tmp = m_mesh2d[meshname]; 01056 } 01057 if(tmp==NULL){ 01058 throw runtime_error("unknown mesh name in TSubLattice<T>::addMesh2DIG:" + meshname); 01059 } 01060 // switch on type,extract params & construc new Mesh2DIG 01061 if(type=="Elastic") 01062 { 01063 AParallelInteractionStorage* new_pis; 01064 ETriMeshIP tmi; 01065 tmi.k=param_buffer.pop_double(); 01066 new_pis = new Mesh2D_PIS_NE<T,EMesh2DInteraction>(tmp,m_ppa,tmi); 01067 m_dpis.insert(make_pair(name,new_pis)); 01068 } else { // unknown type-> throw 01069 throw runtime_error("unknown type in TSubLattice<T>::addMesh2DIG:" + type); 01070 } 01071 01072 01073 console.XDebug() << "end TSubLattice<T>::addTriMeshIG()\n"; 01074 } 01075 01080 template <class T> 01081 void TSubLattice<T>::addBondedMesh2DIG() 01082 { 01083 console.XDebug() << "TSubLattice<T>::addBondedMesh2DIG()\n"; 01084 01085 // MPI buffer 01086 CVarMPIBuffer param_buffer(m_comm); 01087 01088 // receive params 01089 BMesh2DIP param; 01090 param_buffer.receiveBroadcast(0); 01091 01092 // extract name.meshname & params 01093 string name = param_buffer.pop_string(); 01094 console.XDebug() << "BMesh2DIG name: " << name.c_str() << "\n"; 01095 string meshname = param_buffer.pop_string(); 01096 console.XDebug() << "BMesh2DIG mesh name: " << meshname.c_str() << "\n"; 01097 param.k = param_buffer.pop_double(); 01098 console.XDebug() << "BMesh2DIG k : " << param.k << "\n"; 01099 param.brk = param_buffer.pop_double(); 01100 console.XDebug() << "BMesh2DIG r_break: " << param.brk << "\n"; 01101 string buildtype = param_buffer.pop_string(); 01102 console.XDebug() << "BMesh2DIG build type: " << buildtype.c_str() << "\n"; 01103 01104 // get pointer to mesh 01105 Mesh2D* tmp=NULL; 01106 if (m_mesh2d.find(meshname) != m_mesh2d.end()) 01107 { 01108 tmp = m_mesh2d[meshname]; 01109 } 01110 if(tmp==NULL){ 01111 throw runtime_error("unknown mesh name in TSubLattice<T>::addBondedMesh2DIG:" + meshname); 01112 } 01113 01114 // setup new interaction storage 01115 Mesh2D_PIS_EB<T,BMesh2DInteraction>* new_pis=new Mesh2D_PIS_EB<T,BMesh2DInteraction>(tmp,m_ppa,param); 01116 // switch on buildtype, extract buildparam & build 01117 if(buildtype=="BuildByTag"){ 01118 int tag=param_buffer.pop_int(); 01119 int mask=param_buffer.pop_int(); 01120 new_pis->buildFromPPATagged(tag,mask); 01121 m_bpis.insert(make_pair(name,new_pis)); 01122 } else if(buildtype=="BuildByGap"){ 01123 double max_gap=param_buffer.pop_double(); 01124 new_pis->buildFromPPAByGap(max_gap); 01125 m_bpis.insert(make_pair(name,new_pis)); 01126 } else { // unknown type-> throw 01127 throw runtime_error("unknown build type in TSubLattice<T>::addBondedMesh2DIG:" + buildtype); 01128 } 01129 01130 console.XDebug() << "end TSubLattice<T>::addBonded2DMeshIG()\n"; 01131 } 01132 01133 01139 template <class T> 01140 void TSubLattice<T>::addSingleIG() 01141 { 01142 console.XDebug() << "TSubLattice<T>::addSingleIG()\n"; 01143 CVarMPIBuffer param_buffer(m_comm); 01144 01145 // get params 01146 param_buffer.receiveBroadcast(0); 01147 01148 string type = param_buffer.pop_string(); 01149 console.XDebug()<< "SIG type: " << type.c_str() << "\n"; 01150 01151 // setup InteractionGroup 01152 if(type=="Gravity"){ 01153 esys::lsm::BodyForceIGP prms = esys::lsm::BodyForceIGP::extract(¶m_buffer); 01154 01155 // add to map 01156 m_singleParticleInteractions.insert( 01157 std::pair<string,AInteractionGroup<T>*>( 01158 prms.Name(), 01159 new esys::lsm::BodyForceGroup<T>(prms, *m_ppa) 01160 ) 01161 ); 01162 } 01163 else { 01164 throw std::runtime_error( 01165 std::string("Trying to setup SIG of unknown type: ") 01166 + 01167 type 01168 ); 01169 } 01170 01171 console.XDebug() << "end TSubLattice<T>::addSingleIG()\n"; 01172 } 01173 01174 01178 template <class T> 01179 void TSubLattice<T>::addDamping() 01180 { 01181 console.XDebug() << "TSubLattice<T>::addDamping()\n"; 01182 CVarMPIBuffer param_buffer(m_comm); 01183 // get params 01184 param_buffer.receiveBroadcast(0); 01185 01186 string type = param_buffer.pop_string(); 01187 console.XDebug()<< "Damping type: " << type.c_str() << "\n"; 01188 01189 // setup InteractionGroup 01190 doAddDamping(type,param_buffer); 01191 01192 console.XDebug() << "end TSubLattice<T>::addDamping()\n"; 01193 } 01194 01201 template <class T> 01202 bool TSubLattice<T>::doAddDamping(const string& type,CVarMPIBuffer& param_buffer) 01203 { 01204 AParallelInteractionStorage* DG; 01205 string damping_name; 01206 bool res; 01207 01208 if(type=="Damping") 01209 { 01210 CDampingIGP *params=extractDampingIGP(¶m_buffer); 01211 DG=new ParallelInteractionStorage_Single<T,CDamping<T> >(m_ppa,*params); 01212 damping_name="Damping"; 01213 res=true; 01214 } else if (type=="ABCDamping"){ 01215 ABCDampingIGP *params=extractABCDampingIGP(¶m_buffer); 01216 DG=new ParallelInteractionStorage_Single<T,ABCDamping<T> >(m_ppa,*params); 01217 damping_name=params->getName(); 01218 res=true; 01219 } else if (type=="LocalDamping"){ 01220 CLocalDampingIGP *params=extractLocalDampingIGP(¶m_buffer); 01221 DG=new ParallelInteractionStorage_Single<T,CLocalDamping<T> >(m_ppa,*params); 01222 damping_name=params->getName(); 01223 res=true; 01224 }else { 01225 std::stringstream msg; 01226 msg << "Trying to setup Damping of unknown type: " << type; 01227 console.Error() << msg.str() << "\n"; 01228 throw std::runtime_error(msg.str()); 01229 res=false; 01230 } 01231 01232 // add to map 01233 if(res) { 01234 m_damping.insert(make_pair(damping_name,DG)); 01235 m_damping[damping_name]->update(); 01236 } 01237 return res; 01238 } 01239 01244 template <class T> 01245 void TSubLattice<T>::addBondedIG() 01246 { 01247 console.XDebug() << "TSubLattice<T>::addBondedIG()\n"; 01248 CVarMPIBuffer param_buffer(m_comm); 01249 01250 // get params 01251 CBondedIGP param; 01252 param_buffer.receiveBroadcast(0); 01253 param.tag=param_buffer.pop_int(); 01254 string name = param_buffer.pop_string(); 01255 param.k=param_buffer.pop_double(); 01256 param.rbreak=param_buffer.pop_double(); 01257 param.m_scaling=static_cast<bool>(param_buffer.pop_int()); 01258 01259 console.XDebug() 01260 << "Got BondedIG parameters: " << param.tag 01261 << " " << name.c_str() << " " 01262 << param.k << " " << param.rbreak << "\n"; 01263 // setup InteractionGroup 01264 ParallelInteractionStorage_EB<T,CBondedInteraction> *B_PIS = 01265 new ParallelInteractionStorage_EB<T,CBondedInteraction>(m_ppa,param); 01266 01267 // set unbreakable if rbeak<0 01268 if(param.rbreak<0){ 01269 B_PIS->setUnbreakable(true); 01270 console.XDebug() << "set bpis unbreakable\"n"; 01271 } 01272 01273 vector<int> vi(2,-1); 01274 for(size_t i=0;i<m_temp_conn[param.tag].size();i+=2) 01275 { 01276 vi[0] = (m_temp_conn[param.tag][i]); 01277 vi[1] = (m_temp_conn[param.tag][i+1]); 01278 B_PIS->tryInsert(vi); 01279 } 01280 01281 // add InteractionGroup to map 01282 m_bpis.insert(make_pair(name,B_PIS)); 01283 01284 console.XDebug() << "end TSubLattice<T>::addBondedIG()\n"; 01285 } 01286 01291 template <class T> 01292 void TSubLattice<T>::addCappedBondedIG() 01293 { 01294 console.XDebug() << "TSubLattice<T>::addCappedBondedIG()\n"; 01295 CVarMPIBuffer param_buffer(m_comm); 01296 01297 // get params 01298 param_buffer.receiveBroadcast(0); 01299 int tag=param_buffer.pop_int(); 01300 string name = param_buffer.pop_string(); 01301 double k=param_buffer.pop_double(); 01302 double rbreak=param_buffer.pop_double(); 01303 double maxforce=param_buffer.pop_double(); 01304 01305 console.XDebug() 01306 << "Got CappedBondedIG parameters: " << tag 01307 << " " << name.c_str() << " " 01308 << k << " " << rbreak << " " << maxforce << "\n"; 01309 // setup InteractionGroup 01310 CCappedBondedIGP param; 01311 param.k=k; 01312 param.rbreak=rbreak; 01313 param.tag = tag; 01314 param.m_force_limit=maxforce; 01315 ParallelInteractionStorage_EB<T,CCappedBondedInteraction> *B_PIS = 01316 new ParallelInteractionStorage_EB<T,CCappedBondedInteraction>(m_ppa,param); 01317 01318 // set unbreakable if rbeak<0 01319 if(rbreak<0){ 01320 B_PIS->setUnbreakable(true); 01321 console.XDebug() << "set bpis unbreakable\"n"; 01322 } 01323 // recv broadcast connection data 01324 /*console.XDebug() 01325 << "rank=" << m_tml_comm.rank() 01326 << "TSubLattice<T>::addCappedBondedIG(): receiving conn_data.\n"; 01327 01328 vector<int> conn_data; 01329 m_tml_comm.recv_broadcast_cont(conn_data,0); 01330 console.XDebug() 01331 << "rank=" << m_tml_comm.rank() 01332 << "TSubLattice<T>::addBondedIG(): conn_data.size()=" 01333 << conn_data.size() << "\n"; 01334 */ 01335 vector<int> vi(2,-1); 01336 for(size_t i=0;i<m_temp_conn[tag].size();i+=2) 01337 { 01338 vi[0] = (m_temp_conn[tag][i]); 01339 vi[1] = (m_temp_conn[tag][i+1]); 01340 B_PIS->tryInsert(vi); 01341 } 01342 01343 // add InteractionGroup to map 01344 m_bpis.insert(make_pair(name,B_PIS)); 01345 01346 console.XDebug() << "end TSubLattice<T>::addCappedBondedIG()\n"; 01347 } 01348 01349 template <class T> 01350 void TSubLattice<T>::addRotBondedIG() 01351 { 01352 console.Error() << "TSubLattice<T>::addRotBondedIG() => trying to add rotational bonded IG to nonrotational model\n"; 01353 } 01354 01355 template <class T> 01356 void TSubLattice<T>::addRotThermBondedIG() 01357 { 01358 console.Error() << "TSubLattice<T>::addRotThermBondedIG() => trying to add rotational thermal bonded IG to nonrotational model\n"; 01359 } 01360 01365 template <class T> 01366 void TSubLattice<T>::addShortBondedIG() 01367 { 01368 console.XDebug() << "TSubLattice<T>::addShortBondedIG()\n"; 01369 CVarMPIBuffer param_buffer(m_comm); 01370 01371 // get params 01372 param_buffer.receiveBroadcast(0); 01373 int tag=param_buffer.pop_int(); 01374 string name = param_buffer.pop_string(); 01375 double k=param_buffer.pop_double(); 01376 double rbreak=param_buffer.pop_double(); 01377 01378 console.XDebug() 01379 << "Got ShortBondedIG parameters: " << tag 01380 << " " << name.c_str() << " " 01381 << k << " " << rbreak << "\n"; 01382 // setup InteractionGroup 01383 CBondedIGP param; 01384 param.k=k; 01385 param.rbreak=rbreak; 01386 param.tag = tag; 01387 ParallelInteractionStorage_EB<T,CShortBondedInteraction> *B_PIS = 01388 new ParallelInteractionStorage_EB<T,CShortBondedInteraction>(m_ppa,param); 01389 01390 // recv broadcast connection data 01391 /*console.XDebug() 01392 << "rank=" << m_tml_comm.rank() 01393 << "TSubLattice<T>::addShortBondedIG(): receiving conn_data.\n"; 01394 01395 vector<int> conn_data; 01396 m_tml_comm.recv_broadcast_cont(conn_data,0); 01397 console.XDebug() 01398 << "rank=" << m_tml_comm.rank() 01399 << "TSubLattice<T>::addShortBondedIG(): conn_data.size()=" 01400 << conn_data.size() << "\n";*/ 01401 01402 vector<int> vi(2,-1); 01403 for(size_t i=0;i<m_temp_conn[param.tag].size();i+=2) 01404 { 01405 vi[0] = (m_temp_conn[param.tag][i]); 01406 vi[1] = (m_temp_conn[param.tag][i+1]); 01407 B_PIS->tryInsert(vi); 01408 } 01409 01410 // add InteractionGroup to map 01411 m_bpis.insert(make_pair(name,B_PIS)); 01412 01413 console.XDebug() << "end TSubLattice<T>::addShortBondedIG()\n"; 01414 } 01415 01421 template <class T> 01422 void TSubLattice<T>::setExIG() 01423 { 01424 console.XDebug() << "TSubLattice<T>::addExIG()\n"; 01425 CVarMPIBuffer pbuffer(m_comm); 01426 01427 // get params 01428 pbuffer.receiveBroadcast(0); 01429 string s1 = pbuffer.pop_string(); 01430 string s2 = pbuffer.pop_string(); 01431 01432 //console.XDebug()<< s1.c_str() << " " << s2.c_str() << "\n"; 01433 map<string,AParallelInteractionStorage*>::iterator bonded_ig=m_bpis.find(s1); 01434 map<string,AParallelInteractionStorage*>::iterator dynamic_ig=m_dpis.find(s2); 01435 if((bonded_ig!=m_bpis.end())&&(dynamic_ig!=m_dpis.end())) 01436 { 01437 // map iterators point to [key,value] pairs, therefore it->second 01438 // is the pointer to the PIS here 01439 dynamic_ig->second->addExIG(bonded_ig->second); 01440 } 01441 else 01442 { 01443 console.Error() << "TSubLattice<T>::setExIG() - nonexisting interaction group \n"; 01444 } 01445 01446 console.XDebug() << "end TSubLattice<T>::addExIG()\n"; 01447 } 01448 01454 template <class T> 01455 void TSubLattice<T>::removeIG() 01456 { 01457 console.XDebug() << "TSubLattice<T>::removeIG()\n"; 01458 CVarMPIBuffer pbuffer(m_comm); 01459 01460 // get params 01461 pbuffer.receiveBroadcast(0); 01462 string igname = pbuffer.pop_string(); 01463 map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.find(igname); 01464 if(iter!=m_dpis.end()){ 01465 delete m_dpis[igname]; 01466 m_dpis.erase(iter); 01467 } else { 01468 console.Error() << "TSubLattice<T>::removeIG() - nonexisting interaction group - ignore removal\n"; 01469 } 01470 } 01471 01472 01476 template <class T> 01477 void TSubLattice<T>::exchangePos() 01478 { 01479 console.XDebug() << "TSubLattice<T>::exchangePos() \n" ; 01480 01481 m_ppa->exchange(&T::getExchangeValues,&T::setExchangeValues); 01482 01483 console.XDebug() << "end TSubLattice<T>::exchangePos()\n" ; 01484 } 01485 01489 template <class T> 01490 void TSubLattice<T>::zeroForces() 01491 { 01492 console.XDebug() << "TSubLattice<T>::zeroForces()\n"; 01493 01494 // particles 01495 m_ppa->forAllParticles(&T::zeroForce); 01496 // trimeshes 01497 for(map<string,TriMesh*>::iterator iter=m_mesh.begin(); 01498 iter!=m_mesh.end(); 01499 iter++){ 01500 (iter->second)->zeroForces(); 01501 } 01502 01503 // 2d meshes 01504 for(map<string,Mesh2D*>::iterator iter=m_mesh2d.begin(); 01505 iter!=m_mesh2d.end(); 01506 iter++){ 01507 (iter->second)->zeroForces(); 01508 } 01509 // walls 01510 for(typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.begin();iter!=m_WIG.end();iter++) 01511 { 01512 (iter->second)->zeroForce(); 01513 } 01514 console.XDebug() << "end TSubLattice<T>::zeroForces() \n"; 01515 } 01516 01523 template <class T> 01524 void TSubLattice<T>::calcForces() 01525 { 01526 console.XDebug() << "TSubLattice<T>::calcForces() \n"; 01527 01528 // the walls 01529 for(typename map<string,AWallInteractionGroup<T>*>::iterator it=m_WIG.begin();it!=m_WIG.end();it++) 01530 { 01531 (it->second)->calcForces(); 01532 } 01533 // single particle IGs 01534 for( 01535 typename NameIGroupMap::iterator siter=this->m_singleParticleInteractions.begin(); 01536 siter != m_singleParticleInteractions.end(); 01537 siter++ 01538 ) 01539 { 01540 (siter->second)->calcForces(); 01541 } 01542 // dynamically created IGs 01543 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();iter!=m_dpis.end();iter++) 01544 { 01545 (iter->second)->calcForces(); 01546 } 01547 // bonded IGs 01548 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();iter!=m_bpis.end();iter++) 01549 { 01550 (iter->second)->calcForces(); 01551 } 01552 // last, but not least - damping 01553 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();iter!=m_damping.end();iter++) 01554 { 01555 (iter->second)->calcForces(); 01556 } 01557 01558 console.XDebug() << "end TSubLattice<T>::calcForces() \n"; 01559 } 01560 01567 template <class T> 01568 void TSubLattice<T>::setTimeStepSize(double dt) 01569 { 01570 m_dt = dt; 01571 console.XDebug() << "TSubLattice<T>::setTimeStepSize() \n"; 01572 01573 // the walls 01574 for(typename map<string,AWallInteractionGroup<T>*>::iterator it=m_WIG.begin();it!=m_WIG.end();it++) 01575 { 01576 (it->second)->setTimeStepSize(dt); 01577 } 01578 // single particle IGs 01579 for( 01580 typename NameIGroupMap::iterator siter=this->m_singleParticleInteractions.begin(); 01581 siter != m_singleParticleInteractions.end(); 01582 siter++ 01583 ) 01584 { 01585 (siter->second)->setTimeStepSize(dt); 01586 } 01587 // dynamically created IGs 01588 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();iter!=m_dpis.end();iter++) 01589 { 01590 (iter->second)->setTimeStepSize(dt); 01591 } 01592 // bonded IGs 01593 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();iter!=m_bpis.end();iter++) 01594 { 01595 (iter->second)->setTimeStepSize(dt); 01596 } 01597 // last, but not least - damping 01598 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();iter!=m_damping.end();iter++) 01599 { 01600 (iter->second)->setTimeStepSize(dt); 01601 } 01602 01603 console.XDebug() << "end TSubLattice<T>::setTimeStepSize() \n"; 01604 } 01605 01611 template <class T> 01612 void TSubLattice<T>::integrate(double dt) 01613 { 01614 console.XDebug() << "TSubLattice<T>::integrate \n"; 01615 m_ppa->forAllParticles(&T::integrate,dt); 01616 m_ppa->forAllParticles(&T::rescale) ; 01617 console.XDebug() << "end TSubLattice<T>::integrate \n"; 01618 } 01619 01623 template <class T> 01624 void TSubLattice<T>::oneStep() 01625 { 01626 zeroForces(); 01627 calcForces(); 01628 integrate(m_dt); 01629 01630 if (this->getParticleType() == "RotTherm") 01631 { 01632 this->oneStepTherm(); 01633 } 01634 } 01635 01639 template <class T> 01640 void TSubLattice<T>::oneStepTherm() 01641 { 01642 zeroHeat(); // ???? combine? 01643 calcHeatFrict(); 01644 calcHeatTrans(); 01645 integrateTherm(m_dt); 01646 thermExpansion() ; 01647 } 01648 01654 template <class T> 01655 void TSubLattice<T>::integrateTherm(double dt) 01656 { 01657 console.XDebug() << "TSubLattice<T>::integrateTherm \n"; 01658 m_ppa->forAllParticles(&T::integrateTherm,dt); 01659 // m_ppa->forAllParticles(&T::rescale) ; 01660 console.XDebug() << "end TSubLattice<T>::integrateTherm \n"; 01661 } 01662 01663 template <class T> 01664 void TSubLattice<T>::thermExpansion() 01665 { 01666 console.XDebug() << "TSubLattice<T>::thermExpansion() \n"; 01667 m_ppa->forAllParticles(&T::thermExpansion); 01668 // m_ppa->forAllParticles(&T::rescale) ; 01669 console.XDebug() << "end TSubLattice<T>::thermExpansion() \n"; 01670 } 01671 01675 template <class T> 01676 void TSubLattice<T>::zeroHeat() 01677 { 01678 console.XDebug() << "TSubLattice<T>::zeroHeat()\n"; 01679 01680 // particles 01681 m_ppa->forAllParticles(&T::zeroHeat); 01682 01683 /* 01684 // walls 01685 for(typename vector<AWallInteractionGroup<T>*>::iterator iter=m_WIG.begin();iter!=m_WIG.end();iter++) 01686 { 01687 (*iter)->zeroForce(); 01688 } 01689 */ 01690 console.XDebug() << "end TSubLattice<T>::zeroHeat() \n"; 01691 } 01692 01699 template <class T> 01700 void TSubLattice<T>::calcHeatFrict() 01701 { 01702 console.XDebug() << "TSubLattice<T>::calcHeatFrict() \n"; 01703 01704 // dynamically created IGs 01705 for( 01706 typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin(); 01707 iter!=m_dpis.end(); 01708 iter++ 01709 ) 01710 { 01711 (iter->second)->calcHeatFrict(); 01712 } 01713 01714 console.XDebug() << "end TSubLattice<T>::calcHeatFrict() \n"; 01715 } 01716 01717 template <class T> 01718 void TSubLattice<T>::calcHeatTrans() 01719 { 01720 console.XDebug() << "TSubLattice<T>::calcHeatTrans() \n"; 01721 01722 01723 // dynamically created IGs 01724 for( 01725 typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin(); 01726 iter!=m_dpis.end(); 01727 iter++ 01728 ) 01729 { 01730 (iter->second)->calcHeatTrans(); 01731 } 01732 // bonded IGs 01733 for( 01734 typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin(); 01735 iter!=m_bpis.end(); 01736 iter++ 01737 ) 01738 { 01739 (iter->second)->calcHeatTrans(); 01740 } 01741 01742 console.XDebug() << "end TSubLattice<T>::calcHeatTrans() \n"; 01743 } 01744 01748 template <class T> 01749 void TSubLattice<T>::rebuildParticleArray() 01750 { 01751 m_ppa->rebuild(); 01752 } 01753 01757 template <class T> 01758 void TSubLattice<T>::rebuildInteractions() 01759 { 01760 CMPIBarrier barrier(m_worker_comm); 01761 m_pTimers->start("RebuildInteractions"); 01762 m_pTimers->resume("NeighbourSearch"); 01763 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin(); 01764 iter!=m_bpis.end(); 01765 iter++) 01766 { 01767 console.Debug() << "exchg & rebuild BPIS " << iter->first << " at node " << m_rank << "\n"; 01768 (iter->second)->exchange(); 01769 (iter->second)->rebuild(); 01770 } 01771 01772 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin(); 01773 iter!=m_dpis.end(); 01774 iter++) 01775 { 01776 console.Debug() << "exchg & rebuild DPIS " << iter->first << " at node " << m_rank << "\n"; 01777 (iter->second)->exchange(); 01778 m_pTimers->pause("RebuildInteractions"); 01779 m_pTimers->pause("NeighbourSearch"); 01780 barrier.wait("dpis::exchange"); 01781 m_pTimers->resume("RebuildInteractions"); 01782 m_pTimers->resume("NeighbourSearch"); 01783 (iter->second)->rebuild(); 01784 } 01785 resetDisplacements(); 01786 m_pTimers->stop("RebuildInteractions"); 01787 } 01788 01792 template <class T> 01793 void TSubLattice<T>::searchNeighbors() 01794 { 01795 console.Debug() << "CSubLattice<T>::searchNeighbors()\n"; 01796 CMPIBarrier barrier(getWorkerComm()); 01797 m_pTimers->start("NeighbourSearch"); 01798 m_pTimers->start("RebuildParticleArray"); 01799 rebuildParticleArray(); 01800 m_pTimers->stop("RebuildParticleArray"); 01801 m_pTimers->pause("NeighbourSearch"); 01802 barrier.wait("PPA rebuild"); 01803 rebuildInteractions(); 01804 m_pTimers->stop("NeighbourSearch"); 01805 console.Debug() << "end CSubLattice<T>::searchNeighbors()\n"; 01806 } 01807 01813 template <class T> 01814 void TSubLattice<T>::updateInteractions() 01815 { 01816 console.Debug() << "updateInteractions() \n"; 01817 console.Debug() << "m_ppa->getTimeStamp() " << m_ppa->getTimeStamp() << " m_last_ns " << m_last_ns << "\n"; 01818 bool need_update=false; 01819 01820 m_pTimers->start("UpdateBondedInteractions"); 01821 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin(); 01822 iter!=m_bpis.end(); 01823 iter++) 01824 { 01825 bool n=(iter->second)->update(); 01826 need_update=need_update || n; 01827 } 01828 m_pTimers->stop("UpdateBondedInteractions"); 01829 if((m_ppa->getTimeStamp() > m_last_ns) || need_update) 01830 { 01831 m_pTimers->start("UpdateDynamicInteractions"); 01832 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin(); 01833 iter!=m_dpis.end(); 01834 iter++) 01835 { 01836 bool n=(iter->second)->update(); 01837 need_update=need_update || n; 01838 } 01839 m_pTimers->stop("UpdateDynamicInteractions"); 01840 for(typename map<string,AWallInteractionGroup<T>*>::iterator it=m_WIG.begin(); 01841 it!=m_WIG.end(); 01842 it++) 01843 { 01844 (it->second)->Update(m_ppa); 01845 } 01846 for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin(); 01847 iter!=m_damping.end(); 01848 iter++){ 01849 (iter->second)->update(); 01850 } 01851 m_last_ns=m_ppa->getTimeStamp(); 01852 } 01853 01854 console.Debug() << "end TSubLattice<T>::updateInteractions()\n"; 01855 } 01856 01861 template <class T> 01862 void TSubLattice<T>::checkNeighbors() 01863 { 01864 console.Debug() << "TSubLattice<T>::checkNeighbors()\n"; 01865 CMPISGBufferLeaf buffer(m_comm,0,8); 01866 double mdsqr=0; // squared max. displacement 01867 double alpha=0.5*m_alpha; 01868 double srsqr=alpha*alpha; // squared search range 01869 vector<Vec3> displ; // displacements 01870 int res; // result 0/1 01871 01872 // --- particles --- 01873 // get displacement data 01874 m_ppa->forAllParticlesGet(displ,&T::getDisplacement); 01875 01876 // find maximum particle displacement 01877 vector<Vec3>::iterator it=displ.begin(); 01878 while((it!=displ.end())&&(mdsqr<srsqr)) 01879 { 01880 double sqdisp=(*it)*(*it); 01881 mdsqr = ((mdsqr < sqdisp) ? sqdisp : mdsqr); 01882 it++; 01883 //console.XDebug() << "sq. disp: " << sqdisp << "\n"; 01884 } 01885 console.XDebug() << "max squared displacement " << mdsqr << "alpha^2 = " << srsqr << "\n"; 01886 if (mdsqr>srsqr){ 01887 res=1; 01888 } else { 01889 res=0; 01890 } 01891 01892 // --- mesh --- 01893 // only needed if res==0 01894 if(res==0){ 01895 for(map<string,TriMesh*>::iterator iter=m_mesh.begin(); 01896 iter!=m_mesh.end(); 01897 iter++){ 01898 //std::cerr << "check mesh " << iter->first << std::endl; 01899 if(iter->second->hasMovedBy(alpha)){ 01900 res=1; 01901 } 01902 } 01903 } 01904 buffer.append(res); 01905 buffer.send(); 01906 console.Debug() << "end TSubLattice<T>::checkNeighbors()\n"; 01907 } 01908 01909 01913 template <class T> 01914 void TSubLattice<T>::resetDisplacements() 01915 { 01916 console.Debug() << "slave " << m_rank << " resetDisplacements()\n"; 01917 m_ppa->forAllParticles(&T::resetDisplacement); 01918 for(map<string,TriMesh*>::iterator iter=m_mesh.begin(); 01919 iter!=m_mesh.end(); 01920 iter++){ 01921 iter->second->resetCurrentDisplacement(); 01922 } 01923 console.Debug() << "slave " << m_rank << " end resetDisplacements()\n"; 01924 } 01925 01930 template <class T> 01931 void TSubLattice<T>::moveParticleTo() 01932 { 01933 console.Debug() << "TSubLattice<T>::moveParticleTo()\n"; 01934 CVarMPIBuffer buffer(m_comm); 01935 01936 buffer.receiveBroadcast(0); // get data from master 01937 int tag=buffer.pop_int(); 01938 Vec3 mv=buffer.pop_vector(); 01939 m_ppa->forParticleTag(tag,(void (T::*)(Vec3))(&T::moveToRel),mv); 01940 console.Debug() << "end TSubLattice<T>::moveParticleTo()\n"; 01941 } 01942 01947 template <class T> 01948 void TSubLattice<T>::moveTaggedParticlesBy() 01949 { 01950 console.Debug() << "TSubLattice<T>::moveTaggedParticlesBy()\n"; 01951 CVarMPIBuffer buffer(m_comm); 01952 01953 buffer.receiveBroadcast(0); // get data from master 01954 const int tag = buffer.pop_int(); 01955 const Vec3 dx = buffer.pop_vector(); 01956 m_ppa->forParticleTag(tag, (void (T::*)(Vec3))(&T::moveBy),dx); 01957 console.Debug() << "end TSubLattice<T>::moveTaggedParticlesBy()\n"; 01958 } 01959 01960 01961 template <class T> 01962 void TSubLattice<T>::moveSingleParticleTo(int particleId, const Vec3 &posn) 01963 { 01964 m_ppa->forParticle(particleId, (void (T::*)(Vec3))(&T::moveTo), posn); 01965 } 01966 01971 template <class T> 01972 void TSubLattice<T>::moveSingleNode() 01973 { 01974 console.Debug() << "TSubLattice<T>::moveSingleNode()\n"; 01975 CVarMPIBuffer buffer(m_comm); 01976 01977 buffer.receiveBroadcast(0); // get data from master 01978 string name=string(buffer.pop_string()); 01979 int id=buffer.pop_int(); 01980 Vec3 disp=buffer.pop_vector(); 01981 01982 console.XDebug() << "name :" << name << " id : " << id << " disp " << disp << "\n"; 01983 01984 map<string,TriMesh*>::iterator tm=m_mesh.find(name); 01985 if (tm!=m_mesh.end()){ 01986 (tm->second)->moveNode(id,disp); 01987 } else { 01988 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(name); 01989 if(m2d!=m_mesh2d.end()){ 01990 (m2d->second)->moveNode(id,disp); 01991 } 01992 } 01993 console.Debug() << "end TSubLattice<T>::moveSingleNode()\n"; 01994 } 01995 02000 template <class T> 02001 void TSubLattice<T>::moveTaggedNodes() 02002 { 02003 console.Error() << "TSubLattice<T>::moveTaggedNodes() NOT IMPLEMENTED\n"; 02004 throw 02005 std::runtime_error( 02006 "TSubLattice<T>::moveTaggedNodes() NOT IMPLEMENTED\n" 02007 ); 02008 #if 0 02009 CVarMPIBuffer buffer(m_comm); 02010 02011 buffer.receiveBroadcast(0); // get data from master 02012 string name=string(buffer.pop_string()); 02013 int tag=buffer.pop_int(); 02014 Vec3 disp=buffer.pop_vector(); 02015 #endif 02016 } 02017 02024 template <class T> 02025 void TSubLattice<T>::translateMeshBy(const std::string &meshName, const Vec3 &translation) 02026 { 02027 map<string,TriMesh*>::iterator tm=m_mesh.find(meshName); 02028 if (tm != m_mesh.end()){ 02029 (tm->second)->translateBy(translation); 02030 } 02031 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshName); 02032 if(m2d!=m_mesh2d.end()){ 02033 (m2d->second)->translateBy(translation); 02034 } 02035 } 02036 02037 template <class T> 02038 std::pair<double, int> TSubLattice<T>::findParticleNearestTo(const Vec3 &pt) 02039 { 02040 console.Debug() << "TSubLattice<T>::findParticleNearestTo: enter\n"; 02041 const T *pClosest = NULL; 02042 double minDistSqrd = std::numeric_limits<double>::max(); 02043 02044 typename ParticleArray::ParticleIterator it = 02045 m_ppa->getInnerParticleIterator(); 02046 while (it.hasNext()) 02047 { 02048 const T &p = it.next(); 02049 const double distSqrd = (pt - p.getPos()).norm2(); 02050 if (distSqrd < minDistSqrd) 02051 { 02052 minDistSqrd = distSqrd; 02053 pClosest = &p; 02054 } 02055 } 02056 console.Debug() << "TSubLattice<T>::findParticleNearestTo: exit\n"; 02057 return 02058 ( 02059 (pClosest != NULL) 02060 ? 02061 std::make_pair(sqrt(minDistSqrd), pClosest->getID()) 02062 : 02063 std::make_pair(std::numeric_limits<double>::max(), -1) 02064 ); 02065 } 02066 02070 template <class T> 02071 std::pair<int, Vec3> TSubLattice<T>::getParticlePosn(int particleId) 02072 { 02073 const T *particle = NULL; 02074 typename ParticleArray::ParticleIterator it = 02075 m_ppa->getInnerParticleIterator(); 02076 while (it.hasNext()) 02077 { 02078 const T &p = it.next(); 02079 if (p.getID() == particleId) 02080 { 02081 particle = &p; 02082 } 02083 } 02084 if (particle != NULL) 02085 { 02086 return std::make_pair(particleId, particle->getPos()); 02087 } 02088 return std::make_pair(-1,Vec3::ZERO); 02089 } 02090 02094 template <class T> 02095 void TSubLattice<T>::getParticleData(const IdVector &particleIdVector) 02096 { 02097 console.Debug() 02098 << "TSubLattice<T>::getParticleData: enter\n"; 02099 typedef std::set<int> IdSet; 02100 typedef std::vector<T> ParticleVector; 02101 02102 ParticleVector particleVector; 02103 typename ParticleArray::ParticleIterator it = 02104 m_ppa->getInnerParticleIterator(); 02105 if (particleIdVector.size() > 0) 02106 { 02107 IdSet idSet(particleIdVector.begin(), particleIdVector.end()); 02108 console.Debug() 02109 << "TSubLattice<T>::getParticleData: iterating over particles\n"; 02110 while (it.hasNext()) 02111 { 02112 const T &p = it.next(); 02113 if (idSet.find(p.getID()) != idSet.end()) 02114 { 02115 particleVector.push_back(p); 02116 } 02117 } 02118 } 02119 else 02120 { 02121 m_ppa->getAllInnerParticles(particleVector); 02122 } 02123 console.Debug() 02124 << "TSubLattice<T>::getParticleData:" 02125 << " sending particle data of size " << particleVector.size() << "\n"; 02126 m_tml_comm.send_gather_packed(particleVector, 0); 02127 console.Debug() 02128 << "TSubLattice<T>::getParticleData: exit\n"; 02129 } 02130 02134 template <class T> 02135 void TSubLattice<T>::tagParticleNearestTo() 02136 { 02137 console.Debug() << "TSubLattice<T>::tagParticleNearestTo()\n"; 02138 CVarMPIBuffer buffer(m_comm); 02139 02140 buffer.receiveBroadcast(0); // get data from master 02141 int tag=buffer.pop_int(); 02142 int mask=buffer.pop_int(); 02143 Vec3 pos=buffer.pop_vector(); 02144 02145 // warning - this is ugly 02146 T* part_ptr=m_ppa->getParticlePtrByPosition(pos); 02147 if(part_ptr!=NULL){ 02148 int old_tag=part_ptr->getTag(); 02149 int new_tag=(old_tag & (~mask)) | (tag & mask); 02150 part_ptr->setTag(new_tag); 02151 02152 cout << "pos, realpos: " << pos << " " << part_ptr->getPos() << " old tag, new tag " << old_tag << " " << part_ptr->getTag() << endl; 02153 } 02154 console.Debug() << "end TSubLattice<T>::tagParticleNearestTo()\n"; 02155 } 02156 02161 template <class T> 02162 void TSubLattice<T>::setParticleNonDynamic() 02163 { 02164 console.Debug() << "TSubLattice<T>::setParticleNonDynamic()\n"; 02165 CVarMPIBuffer buffer(m_comm); 02166 02167 buffer.receiveBroadcast(0); // get data from master 02168 int tag=buffer.pop_int(); 02169 m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamic)); 02170 console.Debug() << "end TSubLattice<T>::setParticleNonDynamic()\n"; 02171 } 02172 02177 template <class T> 02178 void TSubLattice<T>::setParticleNonRot() 02179 { 02180 console.Debug() << "TSubLattice<T>::setParticleNonRot()\n"; 02181 CVarMPIBuffer buffer(m_comm); 02182 02183 buffer.receiveBroadcast(0); // get data from master 02184 int tag=buffer.pop_int(); 02185 m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamicRot)); 02186 console.Debug() << "end TSubLattice<T>::setParticleNonRot()\n"; 02187 } 02188 02193 template <class T> 02194 void TSubLattice<T>::setParticleNonTrans() 02195 { 02196 console.Debug() << "TSubLattice<T>::setParticleNonTrans()\n"; 02197 CVarMPIBuffer buffer(m_comm); 02198 02199 buffer.receiveBroadcast(0); // get data from master 02200 int tag=buffer.pop_int(); 02201 m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamicLinear)); 02202 console.Debug() << "end TSubLattice<T>::setParticleNonTrans()\n"; 02203 } 02204 02208 template <class T> 02209 void TSubLattice<T>::setTaggedParticleVel() 02210 { 02211 console.Debug() << "TSubLattice<T>::setTaggedParticleVel()\n"; 02212 CVarMPIBuffer buffer(m_comm); 02213 02214 buffer.receiveBroadcast(0); // get data from master 02215 int tag=buffer.pop_int(); 02216 Vec3 v=buffer.pop_vector(); 02217 m_ppa->forParticleTag(tag,(void (T::*)(Vec3))(&T::setVel),v); 02218 console.XDebug() << "end TSubLattice<T>::setTaggedParticleVel()\n"; 02219 } 02220 02224 template <class T> 02225 void TSubLattice<T>::moveWallBy() 02226 { 02227 console.XDebug() << "TSubLattice<T>::moveWallBy()\n"; 02228 CVarMPIBuffer buffer(m_comm); 02229 02230 buffer.receiveBroadcast(0); // get data from master 02231 string wname=buffer.pop_string(); 02232 Vec3 mv=buffer.pop_vector(); 02233 typename map<string,CWall*>::iterator iter=m_walls.find(wname); 02234 if(iter!=m_walls.end()) 02235 { 02236 (iter->second)->moveBy(mv); 02237 } 02238 } 02239 02243 template <class T> 02244 void TSubLattice<T>::setWallNormal() 02245 { 02246 console.XDebug() << "TSubLattice<T>::setWallNormal()\n"; 02247 CVarMPIBuffer buffer(m_comm); 02248 02249 buffer.receiveBroadcast(0); // get data from master 02250 string wname=buffer.pop_string(); 02251 Vec3 wn=buffer.pop_vector(); 02252 typename map<string,CWall*>::iterator iter=m_walls.find(wname); 02253 if(iter!=m_walls.end()) 02254 { 02255 (iter->second)->setNormal(wn); 02256 } 02257 } 02258 02262 template <class T> 02263 void TSubLattice<T>::applyForceToWall() 02264 { 02265 console.XDebug() << "TSubLattice<T>::applyForceToWall()\n"; 02266 CVarMPIBuffer buffer(m_comm); 02267 02268 buffer.receiveBroadcast(0); // get data from master 02269 string wname=buffer.pop_string(); 02270 Vec3 f=buffer.pop_vector(); 02271 typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.find(wname); 02272 if(iter!=m_WIG.end()) 02273 { 02274 (iter->second)->applyForce(f); 02275 } 02276 } 02277 02282 template <class T> 02283 void TSubLattice<T>::setVelocityOfWall() 02284 { 02285 console.XDebug() << "TSubLattice<T>::setVelocityOfWall()\n"; 02286 CVarMPIBuffer buffer(m_comm); 02287 02288 buffer.receiveBroadcast(0); // get data from master 02289 string wname=buffer.pop_string(); 02290 Vec3 v=buffer.pop_vector(); 02291 typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.find(wname); 02292 if(iter!=m_WIG.end()) 02293 { 02294 (iter->second)->setVelocity(v); 02295 } 02296 } 02297 02301 template <class T> 02302 void TSubLattice<T>::setParticleVelocity() 02303 { 02304 console.Debug() << "TSubLattice<T>::setParticleVelocity()\n"; 02305 CVarMPIBuffer buffer(m_comm); 02306 02307 buffer.receiveBroadcast(0); // get data from master 02308 int id=buffer.pop_int(); 02309 Vec3 mv=buffer.pop_vector(); 02310 m_ppa->forParticle(id,(void (T::*)(Vec3))(&T::setVel),mv); 02311 console.XDebug() << "end TSubLattice<T>::setParticleVelocity()\n"; 02312 } 02313 02317 template <class T> 02318 void TSubLattice<T>::setParticleDensity() 02319 { 02320 console.Debug() << "TSubLattice<T>::setParticleDensity()\n"; 02321 CVarMPIBuffer buffer(m_comm); 02322 02323 buffer.receiveBroadcast(0); // get data from master 02324 int tag=buffer.pop_int(); 02325 int mask=buffer.pop_int(); 02326 double rho=buffer.pop_double(); 02327 m_ppa->forParticleTagMask(tag,mask,(void (T::*)(double))(&T::setDensity),rho); 02328 console.XDebug() << "end TSubLattice<T>::setParticleVelocity()\n"; 02329 } 02330 02331 02335 template <class T> 02336 void TSubLattice<T>::sendDataToMaster() 02337 { 02338 console.Debug() << "TSubLattice<T>::sendDataToMaster()\n"; 02339 vector<Vec3> positions; 02340 vector<double> radii; 02341 vector<int> ids; 02342 02343 m_ppa->forAllParticlesGet(positions,(Vec3 (T::*)() const)(&T::getPos)); 02344 m_ppa->forAllParticlesGet(radii,(double (T::*)() const)(&T::getRad)); 02345 m_ppa->forAllParticlesGet(ids,(int (T::*)() const)(&T::getID)); 02346 02347 m_tml_comm.send_gather(positions,0); 02348 m_tml_comm.send_gather(radii,0); 02349 m_tml_comm.send_gather(ids,0); 02350 02351 console.Debug() << "end TSubLattice<T>::sendDataToMaster()\n"; 02352 } 02353 02357 template <class T> 02358 void TSubLattice<T>::countParticles() 02359 { 02360 console.Debug()<<"TSubLattice<T>::countParticles()\n"; 02361 CMPIVarSGBufferLeaf buffer(m_comm,0); 02362 //-- pack particles into buffer 02363 buffer.append(m_ppa->size()); 02364 // send 02365 buffer.send(); 02366 } 02367 02371 template <class T> 02372 void TSubLattice<T>::printStruct() 02373 { 02374 cout<< "My Rank : " << m_rank << "\n" ; 02375 if(m_ppa!=NULL) 02376 { 02377 cout << *m_ppa << endl; 02378 } 02379 } 02380 02381 template <class T> 02382 void TSubLattice<T>::printData() 02383 { 02384 cout << "Data: my rank : " << m_rank << "particles : \n" ; 02385 m_ppa->forAllParticles((void (T::*)())(&T::print)); 02386 } 02387 02388 template <class T> 02389 void TSubLattice<T>::printTimes() 02390 { 02391 console.Debug() << "time spent calculating force : " << forcetime << " sec\n"; 02392 console.Debug() << "time spent communicating : " << commtime << " sec\n"; 02393 console.Debug() << "time spent packing : " << packtime << " sec\n"; 02394 console.Debug() << "time spent unpacking : " << unpacktime << " sec\n"; 02395 } 02396 02397 //-------------- FIELD FUNCTIONS ---------------- 02398 02399 02400 02404 template <class T> 02405 void TSubLattice<T>::addScalarParticleField() 02406 { 02407 // cout << "TSubLattice<T>::addScalarParticleField\n"; 02408 string fieldname; 02409 int id,is_tagged; 02410 02411 m_tml_comm.recv_broadcast_cont(fieldname,0); 02412 //cout << "recvd. fieldname: " << fieldname << "\n"; 02413 m_tml_comm.recv_broadcast(id,0); 02414 //cout << "recvd. id: " << id << "\n"; 02415 m_tml_comm.recv_broadcast(is_tagged,0); 02416 //cout << "recvd. is_tagged: " << is_tagged << "\n"; 02417 02418 typename T::ScalarFieldFunction rdf=T::getScalarFieldFunction(fieldname); 02419 ScalarParticleFieldSlave<T> *new_spfs; 02420 if(is_tagged==0) 02421 { 02422 new_spfs=new ScalarParticleFieldSlave<T>(&m_tml_comm,m_ppa,rdf); 02423 } 02424 else 02425 { 02426 int tag,mask; 02427 m_tml_comm.recv_broadcast(tag,0); 02428 console.XDebug() << "recvd. tag: " << tag << "\n"; 02429 m_tml_comm.recv_broadcast(mask,0); 02430 console.XDebug() << "recvd. mask: " << mask << "\n"; 02431 new_spfs=new ScalarParticleFieldSlaveTagged<T>(&m_tml_comm,m_ppa,rdf,tag,mask); 02432 } 02433 m_field_slaves.insert(make_pair(id,new_spfs)); 02434 } 02435 02439 template <class T> 02440 void TSubLattice<T>::addVectorParticleField() 02441 { 02442 console.XDebug() << "TSubLattice<T>::addVectorParticleField\n"; 02443 string fieldname; 02444 int id,is_tagged; 02445 02446 m_tml_comm.recv_broadcast_cont(fieldname,0); 02447 console.XDebug() << "recvd. fieldname: " << fieldname << "\n"; 02448 m_tml_comm.recv_broadcast(id,0); 02449 console.XDebug() << "recvd. id: " << id << "\n"; 02450 m_tml_comm.recv_broadcast(is_tagged,0); 02451 console.XDebug() << "recvd. is_tagged: " << is_tagged << "\n"; 02452 02453 typename T::VectorFieldFunction rdf=T::getVectorFieldFunction(fieldname); 02454 VectorParticleFieldSlave<T> *new_vpfs; 02455 if(is_tagged==0) 02456 { 02457 new_vpfs=new VectorParticleFieldSlave<T>(&m_tml_comm,m_ppa,rdf); 02458 } 02459 else 02460 { 02461 int tag,mask; 02462 m_tml_comm.recv_broadcast(tag,0); 02463 console.XDebug() << "recvd. tag: " << tag << "\n"; 02464 m_tml_comm.recv_broadcast(mask,0); 02465 console.XDebug() << "recvd. mask: " << mask << "\n"; 02466 new_vpfs=new VectorParticleFieldSlaveTagged<T>(&m_tml_comm,m_ppa,rdf,tag,mask); 02467 } 02468 m_field_slaves.insert(make_pair(id,new_vpfs)); 02469 02470 console.Debug() << "end TSubLattice<T>::addVectorParticleField\n"; 02471 } 02472 02473 02477 template <class T> 02478 void TSubLattice<T>::addScalarInteractionField() 02479 { 02480 console.XDebug() << "TSubLattice<T>::addScalarInteractionField\n"; 02481 string fieldname; 02482 string igname; 02483 string igtype; 02484 int id,is_tagged,tag,mask,is_checked; 02485 02486 m_tml_comm.recv_broadcast_cont(fieldname,0); 02487 console.XDebug() << "recvd. fieldname: " << fieldname << "\n"; 02488 m_tml_comm.recv_broadcast(id,0); 02489 console.XDebug() << "recvd. id: " << id << "\n"; 02490 m_tml_comm.recv_broadcast_cont(igname,0); 02491 console.XDebug() << "recvd. interaction group name: " << igname << "\n"; 02492 m_tml_comm.recv_broadcast_cont(igtype,0); 02493 console.XDebug() << "recvd. interaction group name: " << igtype << "\n"; 02494 m_tml_comm.recv_broadcast(is_tagged,0); 02495 console.XDebug() << "recvd. is_tagged: " << is_tagged << "\n"; 02496 02497 // get interaction group 02498 map<string,AParallelInteractionStorage*>::iterator it=m_dpis.find(igname); 02499 if(is_tagged==1) 02500 { 02501 m_tml_comm.recv_broadcast(tag,0); 02502 m_tml_comm.recv_broadcast(mask,0); 02503 } 02504 m_tml_comm.recv_broadcast(is_checked,0); 02505 console.XDebug() << "recvd. is_checked: " << is_checked << "\n"; 02506 02507 if(it!=m_dpis.end()) 02508 { 02509 console.XDebug() << "found interaction group in dynamic\n"; 02510 AFieldSlave* new_sifs; 02511 new_sifs=(it->second)->generateNewScalarFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask); 02512 m_field_slaves.insert(make_pair(id,new_sifs)); 02513 } 02514 else 02515 { 02516 it=m_bpis.find(igname); 02517 if(it!=m_bpis.end()){ 02518 console.XDebug() << "found interaction group in bonded\n"; 02519 AFieldSlave* new_sifs; 02520 new_sifs=(it->second)->generateNewScalarFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask); 02521 m_field_slaves.insert(make_pair(id,new_sifs)); 02522 } 02523 else // not in dynamic or bonded -> try damping 02524 { 02525 //typename map<string,CDampingGroup<T>*>::iterator it2=m_damping.find(igname); 02526 it=m_damping.find(igname); 02527 if(it!=m_damping.end()) // found it in damping 02528 { 02529 AFieldSlave* new_sifs; 02530 new_sifs=(it->second)->generateNewScalarFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask); 02531 m_field_slaves.insert(make_pair(id,new_sifs)); 02532 } 02533 else // still not found -> unknown name -> error 02534 { 02535 cerr << "ERROR : unknown interaction group name in addScalarInteractionField " << endl; 02536 } 02537 } 02538 } 02539 02540 console.XDebug() << "end TSubLattice<T>::addScalarInteractionField\n"; 02541 } 02542 02546 template <class T> 02547 void TSubLattice<T>::addVectorInteractionField() 02548 { 02549 console.Debug() << "TSubLattice<T>::addVectorInteractionField\n"; 02550 string fieldname; 02551 string igname; 02552 string igtype; 02553 int id,is_tagged,tag,mask,is_checked; 02554 02555 m_tml_comm.recv_broadcast_cont(fieldname,0); 02556 console.XDebug() << "recvd. fieldname: " << fieldname << "\n"; 02557 m_tml_comm.recv_broadcast(id,0); 02558 console.XDebug() << "recvd. id: " << id << "\n"; 02559 m_tml_comm.recv_broadcast_cont(igname,0); 02560 console.XDebug() << "recvd. interaction group name: " << igname << "\n"; 02561 m_tml_comm.recv_broadcast_cont(igtype,0); 02562 console.XDebug() << "recvd. interaction group type: " << igtype << "\n"; 02563 m_tml_comm.recv_broadcast(is_tagged,0); 02564 console.XDebug() << "recvd. is_tagged: " << is_tagged << "\n"; 02565 02566 // get interaction group 02567 map<string,AParallelInteractionStorage*>::iterator it=m_dpis.find(igname); 02568 if(is_tagged==1) 02569 { 02570 m_tml_comm.recv_broadcast(tag,0); 02571 m_tml_comm.recv_broadcast(mask,0); 02572 } 02573 m_tml_comm.recv_broadcast(is_checked,0); 02574 console.XDebug() << "recvd. is_checked: " << is_checked << "\n"; 02575 02576 if(it!=m_dpis.end()) 02577 { 02578 console.XDebug() << "found interaction group in dynamic\n"; 02579 AFieldSlave* new_sifs; 02580 new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask); 02581 if(new_sifs!=NULL){ 02582 m_field_slaves.insert(make_pair(id,new_sifs)); 02583 } else { 02584 console.Error()<<"ERROR: could not generate Field Slave for field " << fieldname << "\n"; 02585 } 02586 } 02587 else 02588 { 02589 it=m_bpis.find(igname); 02590 if(it!=m_bpis.end()){ 02591 console.XDebug() << "found interaction group in bonded\n"; 02592 AFieldSlave* new_sifs; 02593 new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask); 02594 m_field_slaves.insert(make_pair(id,new_sifs)); 02595 } 02596 else // not in dynamic or bonded -> try damping 02597 { 02598 //typename map<string,CDampingGroup<T>*>::iterator it2=m_damping.find(igname); 02599 it=m_damping.find(igname); 02600 if(it!=m_damping.end()) // found it in damping 02601 { 02602 AFieldSlave* new_sifs; 02603 new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask); 02604 m_field_slaves.insert(make_pair(id,new_sifs)); 02605 } 02606 else // still not found -> unknown name -> error 02607 { 02608 cerr << "ERROR : unknown interaction group name in addScalarInteractionField " << endl; 02609 } 02610 } 02611 } 02612 02613 console.Debug() << "end TSubLattice<T>::addVectorInteractionField\n"; 02614 } 02615 02620 template <class T> 02621 void TSubLattice<T>::addVectorTriangleField() 02622 { 02623 console.Debug() << "TSubLattice<T>::addVectorTriangleField()\n"; 02624 string fieldname; 02625 string meshname; 02626 int id; 02627 02628 // receive params 02629 m_tml_comm.recv_broadcast_cont(fieldname,0); 02630 console.XDebug() << "recvd. fieldname: " << fieldname << "\n"; 02631 m_tml_comm.recv_broadcast_cont(meshname,0); 02632 console.XDebug() << "recvd. meshname: " << meshname << "\n"; 02633 m_tml_comm.recv_broadcast(id,0); 02634 console.XDebug() << "recvd. id: " << id << "\n"; 02635 02636 map<string,TriMesh*>::iterator tm=m_mesh.find(meshname); 02637 // if meshname is in trimesh map 02638 if (tm!=m_mesh.end()){ 02639 // get reader function 02640 Triangle::VectorFieldFunction rdf=Triangle::getVectorFieldFunction(fieldname); 02641 // check it 02642 if(rdf==NULL){ 02643 console.Critical() << "NULL rdf !!!\n"; 02644 } 02645 VectorTriangleFieldSlave* new_vfs=new VectorTriangleFieldSlave(&m_tml_comm,tm->second,rdf); 02646 m_field_slaves.insert(make_pair(id,new_vfs)); 02647 } else { 02648 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname); 02649 if(m2d!=m_mesh2d.end()){ 02650 Edge2D::VectorFieldFunction rdf=Edge2D::getVectorFieldFunction(fieldname); 02651 // check it 02652 if(rdf==NULL){ 02653 console.Critical() << "NULL rdf !!!\n"; 02654 } 02655 VectorEdge2DFieldSlave* new_efs= new VectorEdge2DFieldSlave(&m_tml_comm,m2d->second,rdf); 02656 m_field_slaves.insert(make_pair(id,new_efs)); 02657 } 02658 } 02659 console.Debug() << "end TSubLattice<T>::addVectorTriangleField()\n"; 02660 } 02661 02665 template <class T> 02666 void TSubLattice<T>::addScalarTriangleField() 02667 { 02668 console.Debug() << "TSubLattice<T>::addScalarTriangleField()\n"; 02669 string fieldname; 02670 string meshname; 02671 int id; 02672 02673 // receive params 02674 m_tml_comm.recv_broadcast_cont(fieldname,0); 02675 console.XDebug() << "recvd. fieldname: " << fieldname << "\n"; 02676 m_tml_comm.recv_broadcast_cont(meshname,0); 02677 console.XDebug() << "recvd. meshname: " << meshname << "\n"; 02678 m_tml_comm.recv_broadcast(id,0); 02679 console.XDebug() << "recvd. id: " << id << "\n"; 02680 02681 // get reader function 02682 Triangle::ScalarFieldFunction rdf=Triangle::getScalarFieldFunction(fieldname); 02683 // check it 02684 if(rdf==NULL){ 02685 console.Critical() << "NULL rdf !!!\n"; 02686 } 02687 ScalarTriangleFieldSlave* new_vtfs=new ScalarTriangleFieldSlave(&m_tml_comm,m_mesh[meshname],rdf); 02688 m_field_slaves.insert(make_pair(id,new_vtfs)); 02689 console.Debug() << "end TSubLattice<T>::addScalarTriangleField()\n"; 02690 } 02691 02695 template <class T> 02696 void TSubLattice<T>::addVectorWallField() 02697 { 02698 console.XDebug() << "begin TSubLattice<T>::addVectorWallField()\n"; 02699 string fieldname; 02700 string tmp_wallname; 02701 vector<string> wallnames; 02702 int nwall; 02703 int id; 02704 02705 // receive params 02706 m_tml_comm.recv_broadcast_cont(fieldname,0); 02707 console.XDebug() << "recvd. fieldname: " << fieldname << "\n"; 02708 m_tml_comm.recv_broadcast(nwall,0); 02709 console.XDebug() << "recvd. nwall: " << nwall << "\n"; 02710 for(int i=0;i<nwall;i++){ 02711 m_tml_comm.recv_broadcast_cont(tmp_wallname,0); 02712 wallnames.push_back(tmp_wallname); 02713 console.XDebug() << "recvd. wallname: " << tmp_wallname << "\n"; 02714 tmp_wallname.clear(); 02715 } 02716 m_tml_comm.recv_broadcast(id,0); 02717 console.XDebug() << "recvd. id: " << id << "\n"; 02718 02719 // check validity of 1st wall name 02720 map<string,CWall*>::iterator cwalliter=m_walls.find(*(wallnames.begin())); 02721 if(cwalliter==m_walls.end()){ // 1st wall name invalid -> exit 02722 std::stringstream msg; 02723 msg 02724 << "ERROR in addVectorWallField: wallname '" 02725 << *(wallnames.begin()) << " 'invalid. Valid wall names: "; 02726 for (map<string,CWall*>::const_iterator it = m_walls.begin(); it != m_walls.end(); it++) 02727 { 02728 msg << "'" << it->first << "' "; 02729 } 02730 console.Error() << msg.str() << "\n"; 02731 throw std::runtime_error(msg.str()); 02732 } else { // first wall valid -> try to get slave 02733 // get summation flag from wall 02734 int sumflag=(cwalliter->second)->getFieldSummationFlag(fieldname); 02735 // if process 1, send summation flag back to master 02736 if(m_tml_comm.rank()==1){ 02737 m_tml_comm.send(sumflag,0); 02738 } 02739 m_tml_comm.barrier(); 02740 02741 //field slave 02742 AWallFieldSlave* new_fs=(cwalliter->second)->generateVectorFieldSlave(&m_tml_comm,fieldname); 02743 02744 // try to insert other walls 02745 vector<string>::iterator niter=wallnames.begin(); 02746 if(niter!=wallnames.end()) niter++ ; // jump past 1st wall - already got it 02747 while(niter!=wallnames.end()){ 02748 string wname=*niter; 02749 map<string,CWall*>::iterator iter=m_walls.find(wname); 02750 if(iter==m_walls.end()){ // wall name invalid -> exit 02751 std::stringstream msg; 02752 msg 02753 << "ERROR in addVectorWallField: wallname '" 02754 << wname << " 'invalid"; 02755 for (map<string,CWall*>::const_iterator it = m_walls.begin(); it != m_walls.end(); it++) 02756 { 02757 msg << "'" << it->first << "' "; 02758 } 02759 02760 console.Error() << msg.str() << "\n"; 02761 throw std::runtime_error(msg.str()); 02762 } else { 02763 new_fs->addWall(iter->second); 02764 } 02765 niter++; 02766 } 02767 if(new_fs!=NULL){ 02768 m_field_slaves.insert(make_pair(id,new_fs)); 02769 } else { 02770 console.Error() << "ERROR in addVectorWallField: got NULL Slave\n"; 02771 } 02772 } 02773 02774 console.XDebug() << "end TSubLattice<T>::addVectorWallField()\n"; 02775 } 02776 02780 template <class T> 02781 void TSubLattice<T>::sendFieldData() 02782 { 02783 console.Debug() << "TSubLattice<T>::sendFieldData()\n"; 02784 // receive id's of field to send 02785 int id; 02786 m_tml_comm.recv_broadcast(id,0); 02787 console.Debug() << "received field id " << id << " for data collection\n" ; 02788 if(m_field_slaves[id] != NULL) 02789 { 02790 m_field_slaves[id]->sendData(); 02791 } 02792 else 02793 { // can not happen 02794 cerr << "NULL pointer in m_field_slaves!" << endl; 02795 } 02796 // call the sending function of the field 02797 console.Debug() << "end TSubLattice<T>::sendFieldData()\n"; 02798 } 02799 02800 02801 // ---- Checkpointing ---------- 02805 template <class T> 02806 void TSubLattice<T>::saveSnapShotData(std::ostream &oStream) 02807 { 02808 // get precision of output stream and set it to 9 significant digits 02809 std::streamsize oldprec=oStream.precision(9); 02810 02811 // 02812 // Save particle check-point data 02813 // 02814 ParticleArray &particleArray = dynamic_cast<ParticleArray &>(*m_ppa); 02815 typename ParticleArray::ParticleIterator 02816 particleIt(particleArray.getInnerParticleIterator()); 02817 02818 const std::string delim = "\n"; 02819 02820 oStream << particleIt.getNumRemaining() << delim; 02821 while (particleIt.hasNext()) { 02822 particleIt.next().saveSnapShotData(oStream); 02823 oStream << delim; 02824 } 02825 02826 // 02827 // Save Bonded interaction check-point data. 02828 // 02829 typedef std::map<string,AParallelInteractionStorage*> NameBondedInteractionsMap; 02830 typename NameBondedInteractionsMap::iterator it; 02831 oStream << m_bpis.size() << delim; 02832 for (it = m_bpis.begin(); it != m_bpis.end(); it++) { 02833 it->second->saveSnapShotData(oStream); 02834 oStream << delim; 02835 } 02836 02837 // dump trimeshdata (if exists) 02838 oStream << "TMIG " << m_mesh.size() << delim; 02839 for(typename map<string,TriMesh*>::iterator tm_iter=m_mesh.begin(); 02840 tm_iter!=m_mesh.end(); 02841 tm_iter++){ 02842 oStream << tm_iter->first << delim; 02843 tm_iter->second->writeCheckPoint(oStream,delim); 02844 } 02845 02846 // restore output stream to old precision 02847 oStream.precision(oldprec); 02848 } 02849 02853 template <class T> 02854 void TSubLattice<T>::saveCheckPointData(std::ostream &oStream) 02855 { 02856 // get precision of output stream and set it to 12 significant digits 02857 std::streamsize oldprec=oStream.precision(16); 02858 02859 const std::string delim = "\n"; 02860 02861 // 02862 // Save particle check-point data 02863 // 02864 m_ppa->saveCheckPointData(oStream); 02865 02866 // 02867 // Save Bonded interaction check-point data. 02868 // 02869 typedef std::map<string,AParallelInteractionStorage*> NameBondedInteractionsMap; 02870 typename NameBondedInteractionsMap::iterator it; 02871 oStream << m_bpis.size() << delim; 02872 for (it = m_bpis.begin(); it != m_bpis.end(); it++) { 02873 it->second->saveCheckPointData(oStream); 02874 oStream << delim; 02875 } 02876 02877 // 02878 // Save Non-bonded interaction check-point data 02879 // 02880 int count_save=0; 02881 for(std::map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin(); 02882 iter!=m_dpis.end(); 02883 iter++){ 02884 if(iter->second->willSave()) count_save++; 02885 } 02886 oStream << count_save << delim; 02887 for(std::map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin(); 02888 iter!=m_dpis.end(); 02889 iter++){ 02890 if(iter->second->willSave()) { 02891 iter->second->saveCheckPointData(oStream); 02892 oStream << delim; 02893 } 02894 } 02895 02896 // Save walls (name, pos, normal) 02897 oStream << "Walls " << m_walls.size() << delim; 02898 for(map<string,CWall*>::iterator w_iter=m_walls.begin(); 02899 w_iter!=m_walls.end(); 02900 w_iter++){ 02901 oStream << w_iter->first << delim; 02902 w_iter->second->writeCheckPoint(oStream,delim); 02903 } 02904 02905 // dump trimeshdata (if exists) 02906 oStream << "TriMesh " << m_mesh.size() << delim; 02907 for(typename map<string,TriMesh*>::iterator tm_iter=m_mesh.begin(); 02908 tm_iter!=m_mesh.end(); 02909 tm_iter++){ 02910 oStream << tm_iter->first << delim; 02911 tm_iter->second->writeCheckPoint(oStream,delim); 02912 } 02913 // dump 2D mesh data (if exists) 02914 oStream << "Mesh2D " << m_mesh2d.size() << delim; 02915 for(typename map<string,Mesh2D*>::iterator tm_iter=m_mesh2d.begin(); 02916 tm_iter!=m_mesh2d.end(); 02917 tm_iter++){ 02918 oStream << tm_iter->first << delim; 02919 tm_iter->second->writeCheckPoint(oStream,delim); 02920 } 02921 // restore output stream to old precision 02922 oStream.precision(oldprec); 02923 } 02924 02925 template <class T> 02926 void TSubLattice<T>::loadCheckPointData(std::istream &iStream) 02927 { 02928 // get particles 02929 m_ppa->loadCheckPointData(iStream); 02930 02931 // rebuild neighbor table 02932 CMPIBarrier barrier(getWorkerComm()); 02933 m_ppa->rebuild(); 02934 barrier.wait("PPA rebuild"); 02935 02936 //-- get bonds -- 02937 // get nr. of bonded interaction group in the checkpoint file 02938 int nr_bonded_ig; 02939 iStream >> nr_bonded_ig; 02940 02941 // compare with existing bonded particle interaction storage (bpis) 02942 // barf if not equal 02943 if(nr_bonded_ig!=m_bpis.size()){ 02944 std::cerr << "number of bonded interaction groups differ between snapshot and script!" << std::endl; 02945 } else { // numbers fit -> read data 02946 for (map<string,AParallelInteractionStorage*>::iterator it = m_bpis.begin(); 02947 it != m_bpis.end(); 02948 it++) { // for all interaction groups 02949 it->second->loadCheckPointData(iStream); 02950 } 02951 } 02952 //-- get nonbonded interactions -- 02953 // get nr. of non-bonded interaction group in the checkpoint file 02954 int nr_nonbonded_ig; 02955 iStream >> nr_nonbonded_ig; 02956 02957 // compare with existing non-bonded (dynamic) particle interaction storage (dpis) 02958 // barf if not equal 02959 if(nr_nonbonded_ig!=m_dpis.size()){ 02960 std::cerr << "number of dynamic interaction groups differ between snapshot and script!" << std::endl; 02961 } else { // numbers fit -> read data 02962 for (map<string,AParallelInteractionStorage*>::iterator it = m_dpis.begin(); 02963 it != m_dpis.end(); 02964 it++) { // for all interaction groups 02965 it->second->loadCheckPointData(iStream); 02966 } 02967 } 02968 //--- load walls --- 02969 string token; 02970 iStream >> token; 02971 if(token!="Walls") { // found wrong token -> barf 02972 std::cerr << "expected Walls , got " << token << std::endl; 02973 } 02974 // nr. of walls 02975 int nwalls; 02976 iStream >> nwalls; 02977 // read wall names & data 02978 string wname; 02979 for(int i=0;i<nwalls;i++){ 02980 CWall* new_wall = new CWall(); 02981 iStream >> wname; 02982 new_wall->loadCheckPoint(iStream); 02983 m_walls[wname]=new_wall; 02984 } 02985 // --- load meshes -- 02986 int nmesh; 02987 string mname; 02988 // Trimesh (3D) 02989 iStream >> token; 02990 if(token!="TriMesh") { // found wrong token -> barf 02991 std::cerr << "expected TriMesh , got " << token << std::endl; 02992 } 02993 // nr. of meshes 02994 iStream >> nmesh; 02995 // read wall names & data 02996 for(int i=0;i<nmesh;i++){ 02997 TriMesh* new_tm=new TriMesh(); 02998 iStream >> mname; 02999 new_tm->loadCheckPoint(iStream); 03000 m_mesh.insert(make_pair(mname,new_tm)); 03001 } 03002 // Mesh2D 03003 iStream >> token; 03004 if(token!="Mesh2D") { // found wrong token -> barf 03005 std::cerr << "expected Mesh2D , got " << token << std::endl; 03006 } 03007 // nr. of meshes 03008 iStream >> nmesh; 03009 // read wall names & data 03010 for(int i=0;i<nmesh;i++){ 03011 Mesh2D* new_m2d=new Mesh2D(); 03012 iStream >> mname; 03013 new_m2d->loadCheckPoint(iStream); 03014 m_mesh2d.insert(make_pair(mname,new_m2d)); 03015 } 03016 } 03017 03018 // -- mesh data exchange functions -- 03019 03023 template <class T> 03024 void TSubLattice<T>::getMeshNodeRef() 03025 { 03026 console.XDebug() << "TSubLattice<T>::getMeshNodeRef()\n"; 03027 vector<int> ref_vec; 03028 03029 // MPI buffer 03030 CVarMPIBuffer param_buffer(m_comm); 03031 // receive mesh name 03032 param_buffer.receiveBroadcast(0); 03033 string meshname=param_buffer.pop_string(); 03034 console.XDebug() << "Mesh name: " << meshname << "\n"; 03035 03036 // find mesh & collect node ids into array 03037 map<string,TriMesh*>::iterator tm=m_mesh.find(meshname); 03038 if (tm!=m_mesh.end()){ 03039 for(TriMesh::corner_iterator iter=(tm->second)->corners_begin(); 03040 iter!=(tm->second)->corners_end(); 03041 iter++){ 03042 ref_vec.push_back(iter->getID()); 03043 } 03044 } else { 03045 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname); 03046 if(m2d!=m_mesh2d.end()){ 03047 for(Mesh2D::corner_iterator iter=(m2d->second)->corners_begin(); 03048 iter!=(m2d->second)->corners_end(); 03049 iter++){ 03050 ref_vec.push_back(iter->getID()); 03051 } 03052 } else { 03053 console.Critical() << "ERROR - WRONG MESH NAME IN getMeshNodeRef() !! \n"; 03054 } 03055 } 03056 // send back to master 03057 m_tml_comm.send_gather(ref_vec,0); 03058 03059 console.XDebug() << "end TSubLattice<T>::getMeshNodeRef()\n"; 03060 } 03061 03065 template <class T> 03066 void TSubLattice<T>::getMeshFaceRef() 03067 { 03068 console.XDebug() << "TSubLattice<T>::getMeshFaceRef()\n"; 03069 vector<int> ref_vec; 03070 03071 // MPI buffer 03072 CVarMPIBuffer param_buffer(m_comm); 03073 // receive mesh name 03074 param_buffer.receiveBroadcast(0); 03075 string meshname=param_buffer.pop_string(); 03076 console.XDebug() << "Mesh name: " << meshname << "\n"; 03077 03078 // find mesh & collect node ids into array 03079 map<string,TriMesh*>::iterator tm=m_mesh.find(meshname); 03080 if (tm!=m_mesh.end()){ 03081 for(TriMesh::triangle_iterator iter=(tm->second)->triangles_begin(); 03082 iter!=(tm->second)->triangles_end(); 03083 iter++){ 03084 ref_vec.push_back(iter->getID()); 03085 } 03086 } else { 03087 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname); 03088 if(m2d!=m_mesh2d.end()){ 03089 for(Mesh2D::edge_iterator iter=(m2d->second)->edges_begin(); 03090 iter!=(m2d->second)->edges_end(); 03091 iter++){ 03092 ref_vec.push_back(iter->getID()); 03093 } 03094 } else { 03095 console.Critical() << "ERROR - WRONG MESH NAME IN getMeshNodeRef() !! \n"; 03096 } 03097 } 03098 // send back to master 03099 m_tml_comm.send_gather(ref_vec,0); 03100 03101 console.XDebug() << "end TSubLattice<T>::getMeshNodeRef()\n"; 03102 } 03103 03107 template <class T> 03108 void TSubLattice<T>::getMesh2DStress() 03109 { 03110 console.XDebug() << "TSubLattice<T>::getMesh2DStress()\n"; 03111 // receive mesh name 03112 // MPI buffer 03113 CVarMPIBuffer param_buffer(m_comm); 03114 param_buffer.receiveBroadcast(0); 03115 string meshname=param_buffer.pop_string(); 03116 console.XDebug() << "Mesh name: " << meshname << "\n"; 03117 03118 // find mesh & collect data 03119 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname); 03120 if(m2d!=m_mesh2d.end()){ 03121 vector<pair<int,Vec3> > data_vec; 03122 // get data 03123 data_vec=(m2d->second)->forAllEdgesGetIndexed(&Edge2D::getForceDensity); 03124 03125 // send data to master 03126 m_tml_comm.send_gather(data_vec,0); 03127 } else { 03128 console.Critical() << "ERROR - WRONG MESH NAME IN getMesh2DStress() !! \n"; 03129 } 03130 03131 console.XDebug() << "end TSubLattice<T>::getMesh2DStress()\n"; 03132 } 03133 03137 template <class T> 03138 void TSubLattice<T>::getTriMeshForce() 03139 { 03140 console.XDebug() << "TSubLattice<T>::getTriMeshStress(): enter\n"; 03141 // receive mesh name 03142 // MPI buffers 03143 CVarMPIBuffer param_buffer(m_comm); 03144 param_buffer.receiveBroadcast(0); 03145 const std::string meshName = param_buffer.pop_string(); 03146 console.XDebug() << "Mesh name: " << meshName << "\n"; 03147 03148 // find mesh & collect data 03149 map<string,TriMesh*>::iterator m=m_mesh.find(meshName); 03150 if(m != m_mesh.end()){ 03151 vector<pair<int,Vec3> > data_vec; 03152 // get data 03153 data_vec=(m->second)->forAllTrianglesGetIndexed(&Triangle::getForce); 03154 03155 // send data to master 03156 m_tml_comm.send_gather(data_vec,0); 03157 } else { 03158 std::stringstream msg; 03159 msg << "Invalid mesh name: " << meshName << ". No such triangular mesh."; 03160 throw std::runtime_error(msg.str().c_str()); 03161 } 03162 03163 console.XDebug() << "TSubLattice<T>::getTriMeshStress(): exit\n"; 03164 }