ESyS-Particle  4.0.1
SubLattice.hpp
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 &param,
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(&param_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(&param_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(&param_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(&param_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(&param_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(&param_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(&param_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(&param_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 }