ESyS-Particle  4.0.1
SubLattice.hpp
1 
2 // //
3 // Copyright (c) 2003-2011 by The University of Queensland //
4 // Earth Systems Science Computational Centre (ESSCC) //
5 // http://www.uq.edu.au/esscc //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 // -- project includes -
14 
15 #include "Parallel/SubLattice.h"
16 #include "Parallel/MpiInfo.h"
17 #include "Parallel/mpivbuf.h"
18 #include "Parallel/mpisgvbuf.h"
19 #include "Parallel/mpibarrier.h"
20 #include "Parallel/mpia2abuf.h"
21 #include "Model/BondedInteraction.h"
22 #include "Model/CappedBondedInteraction.h"
23 #include "Model/ShortBondedInteraction.h"
24 #include "Model/DampingIGP.h"
25 #include "Model/Damping.h"
26 #include "Model/RotDamping.h"
27 #include "Model/ABCDampingIGP.h"
28 #include "Model/ABCDamping.h"
29 #include "Model/LocalDampingIGP.h"
30 #include "Model/LocalDamping.h"
31 #include "Model/RotLocalDamping.h"
32 #include "Model/FrictionInteraction.h"
33 #include "Model/FractalFriction.h"
34 #include "Model/AdhesiveFriction.h"
35 #include "Model/VWFrictionInteraction.h"
36 #include "Model/RotFricInteraction.h"
37 #include "Model/RotElasticInteraction.h"
38 #include "Model/RotThermFricInteraction.h"
39 #include "Model/RotThermElasticInteraction.h"
40 #include "Model/RotThermBondedInteraction.h"
41 #include "Model/HertzianElasticInteraction.h"
42 #include "Model/HertzianViscoElasticFrictionInteraction.h"
43 #include "Model/HertzianViscoElasticInteraction.h"
44 #include "Model/LinearDashpotInteraction.h"
45 #include "Model/MeshData.h"
46 #include "Model/MeshData2D.h"
47 #include "Model/ETriMeshInteraction.h"
48 #include "Model/BTriMeshInteraction.h"
49 #include "Model/BMesh2DInteraction.h"
50 #include "Model/EMesh2DInteraction.h"
51 
52 // --- parallel storage includes ---
53 #include "ppa/src/pp_array.h"
54 #include "pis/pi_storage_eb.h"
55 #include "pis/pi_storage_ed.h"
56 #include "pis/pi_storage_ed_t.h"
57 #include "pis/pi_storage_ne.h"
58 #include "pis/pi_storage_ne_t.h"
59 #include "pis/pi_storage_single.h"
60 #include "pis/trimesh_pis.h"
61 #include "pis/trimesh_pis_ne.h"
62 #include "pis/trimesh_pis_eb.h"
63 #include "pis/mesh2d_pis_eb.h"
64 #include "pis/mesh2d_pis_ne.h"
65 
66 // --- field includes ---
67 #include "Fields/ScalarParticleFieldSlaveTagged.h"
68 #include "Fields/VectorParticleFieldSlaveTagged.h"
69 #include "Fields/ScalarInteractionFieldSlaveTagged.h"
70 #include "Fields/ScalarParticleFieldSlaveTagged.h"
71 #include "Fields/ScalarInteractionFieldSlaveTagged.h"
72 #include "Fields/CheckedScalarInteractionFieldSlave.h"
73 #include "Fields/CheckedScalarInteractionFieldSlaveTagged.h"
74 #include "Fields/VectorTriangleFieldSlave.h"
75 #include "Fields/ScalarTriangleFieldSlave.h"
76 #include "Fields/VectorEdge2DFieldSlave.h"
77 
78 #include "Model/BodyForceGroup.h"
79 
80 #include <mpi.h>
81 #include <stdlib.h>
82 #include <stdexcept>
83 
84 // -- STL includes --
85 #include <algorithm>
86 #include <stdexcept>
87 #include <boost/limits.hpp>
88 using std::runtime_error;
89 
90 // -- IO includes --
91 #include <iostream>
92 using std::cerr;
93 using std::flush;
94 using std::endl;
96 
97 //----------------------------------------------
98 // TSubLattice functions
99 //----------------------------------------------
100 
108 template <class T>
110  const CLatticeParam &param,
111  int rank,
112  MPI_Comm comm,
113  MPI_Comm worker_comm
114 )
115  : m_ppa(NULL),
116  m_dpis(),
117  m_bpis(),
118  m_singleParticleInteractions(),
119  m_damping(),
120  m_WIG(),
121  m_mesh(),
122  m_mesh2d(),
123  m_dt(0),
124  m_nrange(0),
125  m_alpha(0),
126  m_last_ns(0),
127  m_temp_conn(),
128  m_rank(0),
129  m_comm(MPI_COMM_NULL),
130  m_tml_comm(MPI_COMM_NULL),
131  m_worker_comm(MPI_COMM_NULL),
132  m_tml_worker_comm(MPI_COMM_NULL),
133  m_dims(3, 0),
134  packtime(0),
135  unpacktime(0),
136  commtime(0.0),
137  forcetime(0.0),
138  m_field_slaves(),
139  m_pTimers(NULL)
140 {
141  // cout << "TSubLattice<T>::TSubLattice at " << rank << endl << flush;
142  // -- MPI stuff --
143  m_rank=rank;
144 
145  // set global communicator
146  m_comm=comm;
147  m_tml_comm.setComm(m_comm);
148 
149  m_dims = param.processDims();
150 
151  m_worker_comm=worker_comm;
152  // MPI_Comm_size(m_worker_comm,&m_num_workers);
154 
155 
156  // -- set parameters
157  m_alpha=param.alpha();
158  m_nrange=param.nrange();
159  // cout << "dt,nrange,alpha : " << m_dt << " , " << m_nrange << " , " << m_alpha << "\n";
160 
161  commtime=0.0;
162  packtime=0.0;
163  unpacktime=0.0;
164  forcetime=0.0;
165 
166  m_last_ns=-1;
167 }
168 
169 template <class T>
171 {
172  console.Debug() << "TSubLattice<T>::~TSubLattice(): enter\n";
173  console.Debug()
174  << "TSubLattice<T>::~TSubLattice():"
175  << " deleting wall interaction groups...\n";
176  for(
177  typename map<string,AWallInteractionGroup<T>*>::iterator vit=m_WIG.begin();
178  vit!=m_WIG.end();
179  vit++
180  )
181  {
182  delete vit->second;
183  }
184  console.Debug()
185  << "TSubLattice<T>::~TSubLattice():"
186  << " deleting particle array...\n";
187  if (m_ppa != NULL) delete m_ppa;
188  console.Debug() << "TSubLattice<T>::~TSubLattice(): exit\n";
189 }
190 
197 template <class T>
198 void TSubLattice<T>::initNeighborTable(const Vec3& min,const Vec3& max)
199 {
200  console.XDebug() << "TSubLattice<T>::initNeighborTable(" << min << "," << max << ")\n";
201  // make size fit range
202  double xsize=max.X()-min.X();
203  xsize=m_nrange*ceil(xsize/m_nrange);
204  double ysize=max.Y()-min.Y();
205  ysize=m_nrange*ceil(ysize/m_nrange);
206  double zsize=max.Z()-min.Z();
207  zsize=m_nrange*ceil(zsize/m_nrange);
208  Vec3 grow=Vec3(xsize,ysize,zsize)-(max-min); // size increase
209  Vec3 nmin=min-0.5*grow; // distribute symmetrically
210  Vec3 nmax=max+0.5*grow;
211  console.XDebug() << "range=" << m_nrange << ", new min,max: " << nmin << ", " << nmax << "\n";
212 
213  // construct particle array
214  TML_Comm *ntcomm=new TML_Comm(m_worker_comm);
215  m_ppa=new ParallelParticleArray<T>(ntcomm,m_dims,nmin,nmax,m_nrange,m_alpha);
216  //m_ppa=new ParallelParticleArray<T>(ntcomm,3,nmin,nmax,m_nrange);
217 
218  // makeFields(); // put here to make sure ppa is initialized before makeFields
219 
220  console.XDebug() << "end TSubLattice<T>::initNeighborTable\n";
221 }
222 
223 template <class T>
224 void TSubLattice<T>::do2dCalculations(bool do2d)
225 {
226  T::setDo2dCalculations(do2d);
227 }
228 
229 template <class T>
231 {
232  return m_ppa->getInnerSize();
233 }
234 
242 template <class T>
243 void TSubLattice<T>::initNeighborTable(const Vec3& min,const Vec3& max,const vector<bool>& circ)
244 {
245  console.XDebug() << "TSubLattice<T>::initNeighborTable(" << min << "," << max << ") circ\n";
246  double xsize,ysize,zsize;
247  // if dimension is circular, check if size fits range, otherwise make it fit
248  // x - dim
249  if(circ[0])
250  {
251  xsize=max.X()-min.X();
252  if(fabs(xsize/m_nrange-lrint(xsize/m_nrange))>1e-6)
253  {
254  //console.Critical() << "circular x-dimension doesn't fit range !\n";
255  console.Info() << "Circular x-size incompatible with range, adjusting...\n";
256  m_nrange = xsize/floor(xsize/m_nrange);
257  console.Info() << "New range = " << m_nrange << "\n";
258  }
259  //xsize+=2.0*m_nrange; // padding on the circular ends
260  }
261  else
262  {
263  xsize=max.X()-min.X();
264  xsize=m_nrange*ceil(xsize/m_nrange);
265  }
266  // y - dim
267  if(circ[1])
268  {
269  ysize=max.Y()-min.Y();
270  if(fabs(ysize/m_nrange-lrint(ysize/m_nrange))>1e-6)
271  {
272  console.Critical() << "circular y-dimension doesn't fit range !\n";
273  }
274  ysize+=2.0*m_nrange; // padding on the circular ends
275  }
276  else
277  {
278  ysize=max.Y()-min.Y();
279  ysize=m_nrange*ceil(ysize/m_nrange);
280  }
281  // z - dim
282  if(circ[2])
283  {
284  zsize=max.Z()-min.Z();
285  if(fabs(zsize/m_nrange-lrint(zsize/m_nrange))>1e-6)
286  {
287  console.Critical() << "circular z-dimension doesn't fit range !\n";
288  }
289  zsize+=2.0*m_nrange; // padding on the circular ends
290  }
291  else
292  {
293  zsize=max.Z()-min.Z();
294  zsize=m_nrange*ceil(zsize/m_nrange);
295  }
296  Vec3 grow=Vec3(xsize,ysize,zsize)-(max-min); // size increase
297  Vec3 nmin=min-0.5*grow; // distribute symmetrically
298  Vec3 nmax=max+0.5*grow;
299  console.XDebug() << "range, new min, max: " << m_nrange << " " << nmin << nmax << "\n";
300  // construct particle array
301  TML_Comm *ntcomm=new TML_Comm(m_worker_comm);
302  m_ppa=new ParallelParticleArray<T>(ntcomm,m_dims,circ,nmin,nmax,m_nrange,m_alpha);
303 
304  // makeFields(); // put here to make sure ppa is initialized before makeFields
305 
306  console.XDebug() << "end TSubLattice<T>::initNeighborTable (circ)\n";
307 }
308 
315 template <class T>
317 {
318  console.XDebug() << "TSubLattice<T>::receiveParticles: enter\n";
319 
320  vector<T> recv_buffer;
321  CMPIBarrier barrier(m_comm);
322 
323  m_tml_comm.recv_broadcast_cont_packed(recv_buffer,0);
324  console.XDebug() << "recvd " << recv_buffer.size() << " particles \n";
325  m_ppa->insert(recv_buffer);
326 
327  barrier.wait("TSubLattice<T>::receiveParticles");
328 
329  console.XDebug() << "TSubLattice<T>::receiveParticles: exit\n";
330 }
331 
332 
338 template <class T>
340 {
341  console.XDebug() << "TSubLattice<T>::receiveConnections: enter\n";
342 
343  vector<int> recv_buffer;
344  CMPIBarrier barrier(m_comm);
345 
346  m_tml_comm.recv_broadcast_cont_packed(recv_buffer,0);
347  console.XDebug() << "recvd " << recv_buffer.size() << " connections \n";
348  vector<int>::iterator it;
349  for (it = recv_buffer.begin(); it != recv_buffer.end(); it+=3)
350  {
351  if ( (m_ppa->getParticlePtrByIndex( *(it+1)) == NULL ) ||
352  (m_ppa->getParticlePtrByIndex( *(it+2)) == NULL ) )
353  {
354  continue;
355  }
356  m_temp_conn[*(it)].push_back(*(it+1));
357  m_temp_conn[*(it)].push_back(*(it+2));
358  }
359 
360  barrier.wait("TSubLattice<T>::receiveConnections");
361 
362  console.XDebug() << "TSubLattice<T>::receiveConnections: exit\n";
363 }
364 
365 
369 template <class T>
371 {
372  console.XDebug() << "TSubLattice<T>::addWall: enter\n" ;
373  CVarMPIBuffer param_buffer(m_comm);
374  param_buffer.receiveBroadcast(0);
375 
376  string name=param_buffer.pop_string();
377  Vec3 ipos=param_buffer.pop_vector();
378  Vec3 inorm=param_buffer.pop_vector();
379 
380  m_walls[name]=new CWall(ipos,inorm);
381 
382  console.XDebug() << "TSubLattice<T>::addWall: exit\n" ;
383 }
384 
388 template <class T>
390 {
391  console.XDebug() << "TSubLattice<T>::addElasticWIG: enter\n" ;
392  CVarMPIBuffer param_buffer(m_comm);
393 
394  // receive params from master
395  param_buffer.receiveBroadcast(0);
396 
397  CEWallIGP* wigp=extractEWallIGP(&param_buffer);
398 
399  string wallname=wigp->getWallName();
400  map<string,CWall*>::iterator iter=m_walls.find(wallname);
401  if(iter!=m_walls.end()){
402  AWallInteractionGroup<T>* newCEWIG =
404  &m_tml_worker_comm,
405  m_walls[wallname],
406  wigp
407  );
408  newCEWIG->Update(m_ppa);
409  m_WIG.insert(make_pair(wigp->getName(),newCEWIG));
410  } else {
411  std::stringstream msg;
412  msg << "wall name '" << wallname << "' not found in map of walls";
413  throw std::runtime_error(msg.str().c_str());
414  }
415 
416  delete wigp;
417  console.XDebug() << "TSubLattice<T>::addElasticWIG: exit\n" ;
418 }
419 
420 
424 template <class T>
426 {
427  console.XDebug() << "TSubLattice<T>::addBondedWIG: enter\n" ;
428  CVarMPIBuffer param_buffer(m_comm);
429 
430  // receive params from master
431  param_buffer.receiveBroadcast(0);
432 
433  CBWallIGP* wigp=extractBWallIGP(&param_buffer);
434 
435  string wallname=wigp->getWallName();
436  map<string,CWall*>::iterator iter=m_walls.find(wallname);
437  if(iter!=m_walls.end()){
438  AWallInteractionGroup<T>* newCBWIG=new CBWallInteractionGroup<T>(&m_tml_worker_comm,m_walls[wallname],wigp);
439  newCBWIG->Update(m_ppa);
440  m_WIG.insert(make_pair(wigp->getName(),newCBWIG));
441  } else {
442  console.Error() << "wall name " << wallname << " not found in map of walls\n";
443  }
444 
445  delete wigp;
446  console.XDebug() << "TSubLattice<T>::addBondedWIG: exit\n" ;
447 }
448 
452 template <class T>
454 {
455  console.XDebug() << "TSubLattice<T>::addDirBondedWIG: enter\n" ;
456  CVarMPIBuffer param_buffer(m_comm);
457 
458  // receive params from master
459  param_buffer.receiveBroadcast(0);
460 
461  CSoftBWallIGP* wigp=extractSoftBWallIGP(&param_buffer);
462 
463  string wallname=wigp->getWallName();
464  map<string,CWall*>::iterator iter=m_walls.find(wallname);
465  if(iter!=m_walls.end()){
466  AWallInteractionGroup<T>* newCDWIG=new CSoftBWallInteractionGroup<T>(&m_tml_worker_comm,m_walls[wallname],wigp);
467  newCDWIG->Update(m_ppa);
468  m_WIG.insert(make_pair(wigp->getName(),newCDWIG));
469  } else {
470  console.Error() << "wall name " << wallname << " not found in map of walls\n";
471  }
472 
473  delete wigp;
474  console.XDebug() << "TSubLattice<T>::addDirBondedWIG: exit\n" ;
475 }
476 
480 template <class T>
482 {
483  console.XDebug() << "TSubLattice<T>::getWallPosition: enter\n" ;
484  CVarMPIBuffer param_buffer(m_comm);
485  Vec3 pos;
486 
487  // receive params from master
488  param_buffer.receiveBroadcast(0);
489 
490  std::string wname=param_buffer.pop_string();
491  console.XDebug() << "Wall name: " << wname << "\n";
492 
493  // find wall
494  map<string,CWall*>::iterator iter=m_walls.find(wname);
495  if(iter!=m_walls.end()){
496  pos=(iter->second)->getPos();
497  console.XDebug() << "Wall pos: " << pos << "\n";
498  } else {
499  pos=Vec3(0.0,0.0,0.0);
500  }
501 
502  vector<Vec3> vpos;
503  vpos.push_back(pos);
504  m_tml_comm.send_gather(vpos,0);
505  console.XDebug() << "TSubLattice<T>::getWallPosition: exit\n" ;
506 }
507 
511 template <class T>
513 {
514  console.XDebug() << "TSubLattice<T>::getWallForce: enter\n" ;
515  CVarMPIBuffer param_buffer(m_comm);
516  Vec3 force;
517 
518  // receive params from master
519  param_buffer.receiveBroadcast(0);
520 
521  std::string wname=param_buffer.pop_string();
522  console.XDebug() << "Wall name: " << wname << "\n";
523 
524  // find wall
525  map<string,CWall*>::iterator iter=m_walls.find(wname);
526  if(iter!=m_walls.end()){
527  force=(iter->second)->getForce();
528  console.XDebug() << "Wall force: " << force << "\n";
529  } else {
530  force=Vec3(0.0,0.0,0.0);
531  }
532 
533  vector<Vec3> vforce;
534  vforce.push_back(force);
535  m_tml_comm.send_gather(vforce,0);
536  console.XDebug() << "TSubLattice<T>::getWallForce: exit\n" ;
537 }
538 
542 template <class T>
544 {
545  console.XDebug() << "TSubLattice<T>::addViscWIG: enter\n" ;
546  CVarMPIBuffer param_buffer(m_comm);
547 
548  // receive params from master
549  param_buffer.receiveBroadcast(0);
550 
551  CVWallIGP* wigp=extractVWallIGP(&param_buffer);
552 
553  string wallname=wigp->getWallName();
554  map<string,CWall*>::iterator iter=m_walls.find(wallname);
555  if(iter!=m_walls.end()){
556  AWallInteractionGroup<T>* newCVWIG=new CViscWallIG<T>(&m_tml_worker_comm,m_walls[wallname],wigp);
557  newCVWIG->Update(m_ppa);
558  m_WIG.insert(make_pair(wigp->getName(),newCVWIG));
559  } else {
560  console.Error() << "wall name " << wallname << " not found in map of walls\n";
561  }
562 
563  delete wigp;
564  console.XDebug() << "TSubLattice<T>::addViscWIG: exit\n" ;
565 }
566 
570 template <class T>
572 {
573  console.XDebug() << "TSubLattice<T>::addPairIG()\n";
574  CVarMPIBuffer param_buffer(m_comm,2000);
575 
576  // get params
577  param_buffer.receiveBroadcast(0);
578  string type = param_buffer.pop_string();
579  console.XDebug()<< "PIG type: " << type.c_str() << "\n";
580  string name = param_buffer.pop_string();
581  console.XDebug()<< "PIG name: " << name.c_str() << "\n";
582 
583  doAddPIG(name,type,param_buffer,false);
584 
585  console.XDebug() << "end TSubLattice<T>::addPairIG()\n";
586 }
587 
591 template <class T>
593 {
594  console.XDebug() << "TSubLattice<T>::addTaggedPairIG()\n";
595  CVarMPIBuffer param_buffer(m_comm,2000);
596 
597  // get params
598  param_buffer.receiveBroadcast(0);
599  string type = param_buffer.pop_string();
600  console.XDebug()<< "PIG type: " << type.c_str() << "\n";
601  string name = param_buffer.pop_string();
602  console.XDebug()<< "PIG name: " << name.c_str() << "\n";
603 
604  doAddPIG(name,type,param_buffer,true);
605 
606  console.XDebug() << "end TSubLattice<T>::addTaggedPairIG()\n";
607 }
608 
616 template <class T>
617 bool TSubLattice<T>::doAddPIG(const string& name,const string& type,CVarMPIBuffer& param_buffer, bool tagged)
618 {
619  bool res=false;
620  AParallelInteractionStorage* new_pis = NULL;
621 
622  if(type=="Elastic") {
623  CElasticIGP eigp;
624  eigp.m_k=param_buffer.pop_double();
625  eigp.m_scaling=static_cast<bool>(param_buffer.pop_int());
626  if(tagged){
627  int tag1=param_buffer.pop_int();
628  int mask1=param_buffer.pop_int();
629  int tag2=param_buffer.pop_int();
630  int mask2=param_buffer.pop_int();
631  console.XDebug() << "tag1, mask1, tag2, mask2 "
632  << tag1 << " , " << mask1 << " , "
633  << tag2 << " , " << mask2 << "\n";
634  new_pis=new ParallelInteractionStorage_NE_T<T,CElasticInteraction>(m_ppa,eigp,tag1,mask1,tag2,mask2);
635  }else{
637  }
638  res=true;
639  } else if (type=="Friction") {
640  CFrictionIGP figp;
641  figp.k=param_buffer.pop_double();
642  figp.mu=param_buffer.pop_double();
643  figp.k_s=param_buffer.pop_double();
644  figp.dt=param_buffer.pop_double();
645  figp.m_scaling=static_cast<bool>(param_buffer.pop_int());
646  console.XDebug() << "k,mu,k_s,dt: " << figp.k << " , " << figp.mu << " , "
647  << figp.k_s << " , " << figp.dt << "\n";
649  res=true;
650  } else if (type=="AdhesiveFriction") {
652  figp.k=param_buffer.pop_double();
653  figp.mu=param_buffer.pop_double();
654  figp.k_s=param_buffer.pop_double();
655  figp.dt=param_buffer.pop_double();
656  figp.r_cut=param_buffer.pop_double();
657  console.XDebug()
658  << "k,mu,k_s,dt,r_cut: " << figp.k << " , " << figp.mu << " , "
659  << figp.k_s << " , " << figp.dt << " " << figp.r_cut << "\n";
661  res=true;
662  } else if (type=="FractalFriction") {
663  FractalFrictionIGP figp;
664  figp.k=param_buffer.pop_double();
665  figp.mu_0=param_buffer.pop_double();
666  figp.k_s=param_buffer.pop_double();
667  figp.dt=param_buffer.pop_double();
668  console.XDebug() << "k,mu_0,k_s,dt: " << figp.k << " , " << figp.mu_0 << " , "
669  << figp.k_s << " , " << figp.dt << "\n";
670  figp.x0=param_buffer.pop_double();
671  figp.y0=param_buffer.pop_double();
672  figp.dx=param_buffer.pop_double();
673  figp.dy=param_buffer.pop_double();
674  figp.nx=param_buffer.pop_int();
675  figp.ny=param_buffer.pop_int();
676  console.XDebug()
677  <<"x0,y0,dx,dy,nx,ny: "
678  << figp.x0 << " , " << figp.y0 << " , "
679  << figp.dx << " , " << figp.dy << " ,"
680  << figp.nx << " , " << figp.ny << "\n";
681  figp.mu = boost::shared_ptr<double>(new double[figp.nx*figp.ny]);
682 
683  for(int i=0;i<figp.nx*figp.ny;i++)
684  {
685  (figp.mu.get())[i]=param_buffer.pop_double();
686  // console.XDebug() << i << " , " << figp.mu[i] << "\n";
687  }
688  new_pis = new ParallelInteractionStorage_ED<T,CFractalFriction>(m_ppa,figp);
689  res=true;
690  } else if(type=="VWFriction") {
691  VWFrictionIGP figp;
692 
693  figp.k=param_buffer.pop_double();
694  figp.mu=param_buffer.pop_double();
695  figp.k_s=param_buffer.pop_double();
696  figp.dt=param_buffer.pop_double();
697  figp.m_alpha=param_buffer.pop_double();
698  console.XDebug()
699  << "k,mu,k_s,dt,alpha: " << figp.k << " , " << figp.mu << " , "
700  << figp.k_s << " , " << figp.dt << "\n";
701  new_pis=new ParallelInteractionStorage_ED<T,CVWFriction>(m_ppa,figp);
702  res=true;
703  } else if(type=="RotElastic"){
704  CRotElasticIGP reigp;
705  reigp.m_kr=param_buffer.pop_double();
707  res=true;
708  } else if (type=="RotFriction"){
709  CRotFrictionIGP rfigp;
710  rfigp.k=param_buffer.pop_double();
711  rfigp.mu_s=param_buffer.pop_double();
712  rfigp.mu_d=param_buffer.pop_double();
713  rfigp.k_s=param_buffer.pop_double();
714  rfigp.dt=param_buffer.pop_double();
715  rfigp.scaling=static_cast<bool>(param_buffer.pop_int());
716  console.XDebug()
717  << "k,mu_s,mu_d,k_s,dt,scaling: " << rfigp.k << " , "
718  << rfigp.mu_s << " , " << rfigp.mu_d << " , "
719  << rfigp.k_s << " , " << rfigp.dt << " , " << rfigp.scaling << "\n";
720  if(tagged){
721  int tag1=param_buffer.pop_int();
722  int mask1=param_buffer.pop_int();
723  int tag2=param_buffer.pop_int();
724  int mask2=param_buffer.pop_int();
725  console.XDebug() << "tag1, mask1, tag2, mask2 "
726  << tag1 << " , " << mask1 << " , "
727  << tag2 << " , " << mask2 << "\n";
728  new_pis=new ParallelInteractionStorage_ED_T<CRotParticle,CRotFrictionInteraction>(m_ppa,rfigp,tag1,mask1,tag2,mask2);
729  } else {
731  }
732  res=true;
733  } else if (type == "RotThermElastic") {
734  CRotThermElasticIGP eigp;
735  eigp.m_kr = param_buffer.pop_double();
736  eigp.diffusivity = param_buffer.pop_double();
737  console.XDebug()
738  << "k=" << eigp.m_kr << " , "
739  << "diffusivity=" << eigp.diffusivity << "\n";
740 
741  new_pis =
745  >(
746  m_ppa,
747  eigp
748  );
749  res=true;
750  } else if (type == "RotThermFriction") {
751  CRotThermFrictionIGP rfigp;
752  rfigp.k=param_buffer.pop_double();
753  rfigp.mu_s=param_buffer.pop_double();
754  rfigp.mu_d=param_buffer.pop_double();
755  rfigp.k_s=param_buffer.pop_double();
756  rfigp.diffusivity=param_buffer.pop_double();
757  rfigp.dt=param_buffer.pop_double();
758  console.XDebug()
759  << "k=" << rfigp.k << " , "
760  << "mu_d=" << rfigp.mu_d << " , "
761  << "mu_s=" << rfigp.mu_s << " , "
762  << "k_s=" << rfigp.k_s << " , "
763  << "diffusivity=" << rfigp.diffusivity << " , "
764  << "dt=" << rfigp.dt << "\n";
765 
766  new_pis =
770  >(
771  m_ppa,
772  rfigp
773  );
774  res=true;
775  } else if(type=="HertzianElastic") {
776  CHertzianElasticIGP heigp;
777  heigp.m_E=param_buffer.pop_double();
778  heigp.m_nu=param_buffer.pop_double();
779  if(tagged){
780  int tag1=param_buffer.pop_int();
781  int mask1=param_buffer.pop_int();
782  int tag2=param_buffer.pop_int();
783  int mask2=param_buffer.pop_int();
784  new_pis=new ParallelInteractionStorage_NE_T<T,CHertzianElasticInteraction>(m_ppa,heigp,tag1,mask1,tag2,mask2);
785  }else{
787  }
788  res=true;
789  } else if(type=="HertzianViscoElasticFriction") {
791  hvefigp.m_A=param_buffer.pop_double();
792  hvefigp.m_E=param_buffer.pop_double();
793  hvefigp.m_nu=param_buffer.pop_double();
794  hvefigp.mu=param_buffer.pop_double();
795  hvefigp.k_s=param_buffer.pop_double();
796  hvefigp.dt=param_buffer.pop_double();
797  if(tagged){
798  int tag1=param_buffer.pop_int();
799  int mask1=param_buffer.pop_int();
800  int tag2=param_buffer.pop_int();
801  int mask2=param_buffer.pop_int();
802  new_pis=new ParallelInteractionStorage_NE_T<T,CHertzianViscoElasticFrictionInteraction>(m_ppa,hvefigp,tag1,mask1,tag2,mask2);
803  }else{
805  }
806  res=true;
807  } else if(type=="HertzianViscoElastic") {
809  hveigp.m_A=param_buffer.pop_double();
810  hveigp.m_E=param_buffer.pop_double();
811  hveigp.m_nu=param_buffer.pop_double();
812  if(tagged){
813  int tag1=param_buffer.pop_int();
814  int mask1=param_buffer.pop_int();
815  int tag2=param_buffer.pop_int();
816  int mask2=param_buffer.pop_int();
817  new_pis=new ParallelInteractionStorage_NE_T<T,CHertzianViscoElasticInteraction>(m_ppa,hveigp,tag1,mask1,tag2,mask2);
818  }else{
820  }
821  res=true;
822  } else if(type=="LinearDashpot") {
823  CLinearDashpotIGP heigp;
824  heigp.m_damp=param_buffer.pop_double();
825  heigp.m_cutoff=param_buffer.pop_double();
826  if(tagged){
827  int tag1=param_buffer.pop_int();
828  int mask1=param_buffer.pop_int();
829  int tag2=param_buffer.pop_int();
830  int mask2=param_buffer.pop_int();
831  new_pis=new ParallelInteractionStorage_NE_T<T,CLinearDashpotInteraction>(m_ppa,heigp,tag1,mask1,tag2,mask2);
832  }else{
834  }
835  res=true;
836  } else {
837  cerr << "Unknown interaction group name "
838  << type
839  << " in TSubLattice<T>::addPairIG()" << endl;
840  }
841 
842  // add InteractionGroup to map
843  if(res) m_dpis.insert(make_pair(name,new_pis));
844 
845  return res;
846 }
847 
848 
852 template <class T>
854 {
855  console.XDebug() << "TSubLattice<T>::addTriMesh()\n";
856 
857  // MPI buffer
858  CVarMPIBuffer param_buffer(m_comm);
859  // data buffers
860  vector<MeshNodeData> node_recv_buffer;
861  vector<MeshTriData> tri_recv_buffer;
862 
863  // receive params
864  param_buffer.receiveBroadcast(0);
865 
866  // extract name
867  string name = param_buffer.pop_string();
868  console.XDebug()<< "TriMesh name: " << name.c_str() << "\n";
869 
870  // receive nodes
871  m_tml_comm.recv_broadcast_cont_packed(node_recv_buffer,0);
872  console.XDebug() << "recvd " << node_recv_buffer.size() << " nodes \n";
873 
874  // receive triangles
875  m_tml_comm.recv_broadcast_cont_packed(tri_recv_buffer,0);
876  console.XDebug() << "recvd " << tri_recv_buffer.size() << " triangles \n";
877 
878  // load mesh into new TriMesh
879  TriMesh* new_tm=new TriMesh();
880  new_tm->LoadMesh(node_recv_buffer,tri_recv_buffer);
881 
882  m_mesh.insert(make_pair(name,new_tm));
883 
884  console.XDebug() << "end TSubLattice<T>::addTriMesh()\n";
885 }
886 
890 template <class T>
892 {
893  console.XDebug() << "TSubLattice<T>::addTriMeshIG()\n";
894 
895  // MPI buffer
896  CVarMPIBuffer param_buffer(m_comm);
897 
898  // receive params
899  param_buffer.receiveBroadcast(0);
900 
901  // extract name & type
902  string type = param_buffer.pop_string();
903  console.XDebug()<< "TriMeshIG type: " << type.c_str() << "\n";
904  string name = param_buffer.pop_string();
905  console.XDebug()<< "TriMeshIG name: " << name.c_str() << "\n";
906  string meshname = param_buffer.pop_string();
907  console.XDebug()<< "TriMeshIG mesh name: " << meshname.c_str() << "\n";
908 
909  // get pointer to mesh
910  TriMesh* tmp=NULL;
911  if (m_mesh.find(meshname) != m_mesh.end())
912  {
913  tmp = m_mesh[meshname];
914  }
915  if(tmp==NULL){
916  throw runtime_error("unknown mesh name in TSubLattice<T>::addTriMeshIG:" + meshname);
917  }
918  // switch on type,extract params & construc new TriMeshIG
919  if(type=="Elastic")
920  {
922  ETriMeshIP tmi;
923  tmi.k=param_buffer.pop_double();
924  new_pis = new TriMesh_PIS_NE<T,ETriMeshInteraction>(tmp,m_ppa,tmi);
925  m_dpis.insert(make_pair(name,new_pis));
926  } else { // unknown type-> throw
927  throw runtime_error("unknown type in TSubLattice<T>::addTriMeshIG:" + type);
928  }
929 
930 
931  console.XDebug() << "end TSubLattice<T>::addTriMeshIG()\n";
932 }
933 
937 template <class T>
939 {
940  console.XDebug() << "TSubLattice<T>::addBondedTriMeshIG()\n";
941 
942  // MPI buffer
943  CVarMPIBuffer param_buffer(m_comm);
944 
945  // receive params
946  BTriMeshIP param;
947  param_buffer.receiveBroadcast(0);
948 
949  // extract name.meshname & params
950  string name = param_buffer.pop_string();
951  console.XDebug()<< "BTriMeshIG name: " << name.c_str() << "\n";
952  string meshname = param_buffer.pop_string();
953  console.XDebug()<< "BTriMeshIG mesh name: " << meshname.c_str() << "\n";
954  param.k = param_buffer.pop_double();
955  console.XDebug()<< "BTriMeshIG k : " << param.k << "\n";
956  param.brk = param_buffer.pop_double();
957  console.XDebug()<< "BTriMeshIG r_break: " << param.brk << "\n";
958  string buildtype = param_buffer.pop_string();
959  console.XDebug()<< "BTriMeshIG build type: " << buildtype.c_str() << "\n";
960 
961  // get pointer to mesh
962  TriMesh* tmp=NULL;
963  if (m_mesh.find(meshname) != m_mesh.end())
964  {
965  tmp = m_mesh[meshname];
966  }
967  if(tmp==NULL){
968  throw runtime_error("unknown mesh name in TSubLattice<T>::addTriMeshIG:" + meshname);
969  }
970 
971  // setup new interaction storage
973  // switch on buildtype, extract buildparam & build
974  if(buildtype=="BuildByTag"){
975  int tag=param_buffer.pop_int();
976  int mask=param_buffer.pop_int();
977  new_pis->buildFromPPATagged(tag,mask);
978  m_bpis.insert(make_pair(name,new_pis));
979  } else if(buildtype=="BuildByGap"){
980  double max_gap=param_buffer.pop_double();
981  new_pis->buildFromPPAByGap(max_gap);
982  m_bpis.insert(make_pair(name,new_pis));
983  } else { // unknown type-> throw
984  throw runtime_error("unknown build type in TSubLattice<T>::addBondedTriMeshIG:" + buildtype);
985  }
986 
987  console.XDebug() << "end TSubLattice<T>::addBondedTriMeshIG()\n";
988 }
989 
993 template <class T>
995 {
996  console.XDebug() << "TSubLattice<T>::addMesh2D()\n";
997 
998  // MPI buffer
999  CVarMPIBuffer param_buffer(m_comm);
1000  // data buffers
1001  vector<MeshNodeData2D> node_recv_buffer;
1002  vector<MeshEdgeData2D> edge_recv_buffer;
1003 
1004  // receive params
1005  param_buffer.receiveBroadcast(0);
1006 
1007  // extract name
1008  string name = param_buffer.pop_string();
1009  console.XDebug()<< "Mesh2D name: " << name.c_str() << "\n";
1010 
1011  // receive nodes
1012  m_tml_comm.recv_broadcast_cont_packed(node_recv_buffer,0);
1013  console.XDebug() << "recvd " << node_recv_buffer.size() << " nodes \n";
1014 
1015  // receive edges
1016  m_tml_comm.recv_broadcast_cont_packed(edge_recv_buffer,0);
1017  console.XDebug() << "recvd " << edge_recv_buffer.size() << " edges \n";
1018 
1019  // load mesh into new 2D Mesh
1020  Mesh2D* new_tm=new Mesh2D();
1021  new_tm->LoadMesh(node_recv_buffer,edge_recv_buffer);
1022 
1023  m_mesh2d.insert(make_pair(name,new_tm));
1024 
1025  console.XDebug() << "end TSubLattice<T>::addMesh2D()\n";
1026 }
1027 
1032 template <class T>
1034 {
1035 console.XDebug() << "TSubLattice<T>::addMesh2DIG()\n";
1036 
1037  // MPI buffer
1038  CVarMPIBuffer param_buffer(m_comm);
1039 
1040  // receive params
1041  param_buffer.receiveBroadcast(0);
1042 
1043  // extract name & type
1044  string type = param_buffer.pop_string();
1045  console.XDebug()<< "Mesh2DIG type: " << type.c_str() << "\n";
1046  string name = param_buffer.pop_string();
1047  console.XDebug()<< "Mesh2DIG name: " << name.c_str() << "\n";
1048  string meshname = param_buffer.pop_string();
1049  console.XDebug()<< "Mesh2DIG mesh name: " << meshname.c_str() << "\n";
1050 
1051  // get pointer to mesh
1052  Mesh2D* tmp=NULL;
1053  if (m_mesh2d.find(meshname) != m_mesh2d.end())
1054  {
1055  tmp = m_mesh2d[meshname];
1056  }
1057  if(tmp==NULL){
1058  throw runtime_error("unknown mesh name in TSubLattice<T>::addMesh2DIG:" + meshname);
1059  }
1060  // switch on type,extract params & construc new Mesh2DIG
1061  if(type=="Elastic")
1062  {
1063  AParallelInteractionStorage* new_pis;
1064  ETriMeshIP tmi;
1065  tmi.k=param_buffer.pop_double();
1066  new_pis = new Mesh2D_PIS_NE<T,EMesh2DInteraction>(tmp,m_ppa,tmi);
1067  m_dpis.insert(make_pair(name,new_pis));
1068  } else { // unknown type-> throw
1069  throw runtime_error("unknown type in TSubLattice<T>::addMesh2DIG:" + type);
1070  }
1071 
1072 
1073  console.XDebug() << "end TSubLattice<T>::addTriMeshIG()\n";
1074 }
1075 
1080 template <class T>
1082 {
1083  console.XDebug() << "TSubLattice<T>::addBondedMesh2DIG()\n";
1084 
1085  // MPI buffer
1086  CVarMPIBuffer param_buffer(m_comm);
1087 
1088  // receive params
1089  BMesh2DIP param;
1090  param_buffer.receiveBroadcast(0);
1091 
1092  // extract name.meshname & params
1093  string name = param_buffer.pop_string();
1094  console.XDebug() << "BMesh2DIG name: " << name.c_str() << "\n";
1095  string meshname = param_buffer.pop_string();
1096  console.XDebug() << "BMesh2DIG mesh name: " << meshname.c_str() << "\n";
1097  param.k = param_buffer.pop_double();
1098  console.XDebug() << "BMesh2DIG k : " << param.k << "\n";
1099  param.brk = param_buffer.pop_double();
1100  console.XDebug() << "BMesh2DIG r_break: " << param.brk << "\n";
1101  string buildtype = param_buffer.pop_string();
1102  console.XDebug() << "BMesh2DIG build type: " << buildtype.c_str() << "\n";
1103 
1104  // get pointer to mesh
1105  Mesh2D* tmp=NULL;
1106  if (m_mesh2d.find(meshname) != m_mesh2d.end())
1107  {
1108  tmp = m_mesh2d[meshname];
1109  }
1110  if(tmp==NULL){
1111  throw runtime_error("unknown mesh name in TSubLattice<T>::addBondedMesh2DIG:" + meshname);
1112  }
1113 
1114  // setup new interaction storage
1116  // switch on buildtype, extract buildparam & build
1117  if(buildtype=="BuildByTag"){
1118  int tag=param_buffer.pop_int();
1119  int mask=param_buffer.pop_int();
1120  new_pis->buildFromPPATagged(tag,mask);
1121  m_bpis.insert(make_pair(name,new_pis));
1122  } else if(buildtype=="BuildByGap"){
1123  double max_gap=param_buffer.pop_double();
1124  new_pis->buildFromPPAByGap(max_gap);
1125  m_bpis.insert(make_pair(name,new_pis));
1126  } else { // unknown type-> throw
1127  throw runtime_error("unknown build type in TSubLattice<T>::addBondedMesh2DIG:" + buildtype);
1128  }
1129 
1130  console.XDebug() << "end TSubLattice<T>::addBonded2DMeshIG()\n";
1131 }
1132 
1133 
1139 template <class T>
1141 {
1142  console.XDebug() << "TSubLattice<T>::addSingleIG()\n";
1143  CVarMPIBuffer param_buffer(m_comm);
1144 
1145  // get params
1146  param_buffer.receiveBroadcast(0);
1147 
1148  string type = param_buffer.pop_string();
1149  console.XDebug()<< "SIG type: " << type.c_str() << "\n";
1150 
1151  // setup InteractionGroup
1152  if(type=="Gravity"){
1153  esys::lsm::BodyForceIGP prms = esys::lsm::BodyForceIGP::extract(&param_buffer);
1154 
1155  // add to map
1156  m_singleParticleInteractions.insert(
1157  std::pair<string,AInteractionGroup<T>*>(
1158  prms.Name(),
1159  new esys::lsm::BodyForceGroup<T>(prms, *m_ppa)
1160  )
1161  );
1162  }
1163  else {
1164  throw std::runtime_error(
1165  std::string("Trying to setup SIG of unknown type: ")
1166  +
1167  type
1168  );
1169  }
1170 
1171  console.XDebug() << "end TSubLattice<T>::addSingleIG()\n";
1172 }
1173 
1174 
1178 template <class T>
1180 {
1181  console.XDebug() << "TSubLattice<T>::addDamping()\n";
1182  CVarMPIBuffer param_buffer(m_comm);
1183  // get params
1184  param_buffer.receiveBroadcast(0);
1185 
1186  string type = param_buffer.pop_string();
1187  console.XDebug()<< "Damping type: " << type.c_str() << "\n";
1188 
1189  // setup InteractionGroup
1190  doAddDamping(type,param_buffer);
1191 
1192  console.XDebug() << "end TSubLattice<T>::addDamping()\n";
1193 }
1194 
1201 template <class T>
1202 bool TSubLattice<T>::doAddDamping(const string& type,CVarMPIBuffer& param_buffer)
1203 {
1205  string damping_name;
1206  bool res;
1207 
1208  if(type=="Damping")
1209  {
1210  CDampingIGP *params=extractDampingIGP(&param_buffer);
1211  DG=new ParallelInteractionStorage_Single<T,CDamping<T> >(m_ppa,*params);
1212  damping_name="Damping";
1213  res=true;
1214  } else if (type=="ABCDamping"){
1215  ABCDampingIGP *params=extractABCDampingIGP(&param_buffer);
1216  DG=new ParallelInteractionStorage_Single<T,ABCDamping<T> >(m_ppa,*params);
1217  damping_name=params->getName();
1218  res=true;
1219  } else if (type=="LocalDamping"){
1220  CLocalDampingIGP *params=extractLocalDampingIGP(&param_buffer);
1222  damping_name=params->getName();
1223  res=true;
1224  }else {
1225  std::stringstream msg;
1226  msg << "Trying to setup Damping of unknown type: " << type;
1227  console.Error() << msg.str() << "\n";
1228  throw std::runtime_error(msg.str());
1229  res=false;
1230  }
1231 
1232  // add to map
1233  if(res) {
1234  m_damping.insert(make_pair(damping_name,DG));
1235  m_damping[damping_name]->update();
1236  }
1237  return res;
1238 }
1239 
1244 template <class T>
1246 {
1247  console.XDebug() << "TSubLattice<T>::addBondedIG()\n";
1248  CVarMPIBuffer param_buffer(m_comm);
1249 
1250  // get params
1251  CBondedIGP param;
1252  param_buffer.receiveBroadcast(0);
1253  param.tag=param_buffer.pop_int();
1254  string name = param_buffer.pop_string();
1255  param.k=param_buffer.pop_double();
1256  param.rbreak=param_buffer.pop_double();
1257  param.m_scaling=static_cast<bool>(param_buffer.pop_int());
1258 
1259  console.XDebug()
1260  << "Got BondedIG parameters: " << param.tag
1261  << " " << name.c_str() << " "
1262  << param.k << " " << param.rbreak << "\n";
1263  // setup InteractionGroup
1266 
1267  // set unbreakable if rbeak<0
1268  if(param.rbreak<0){
1269  B_PIS->setUnbreakable(true);
1270  console.XDebug() << "set bpis unbreakable\"n";
1271  }
1272 
1273  vector<int> vi(2,-1);
1274  for(size_t i=0;i<m_temp_conn[param.tag].size();i+=2)
1275  {
1276  vi[0] = (m_temp_conn[param.tag][i]);
1277  vi[1] = (m_temp_conn[param.tag][i+1]);
1278  B_PIS->tryInsert(vi);
1279  }
1280 
1281  // add InteractionGroup to map
1282  m_bpis.insert(make_pair(name,B_PIS));
1283 
1284  console.XDebug() << "end TSubLattice<T>::addBondedIG()\n";
1285 }
1286 
1291 template <class T>
1293 {
1294  console.XDebug() << "TSubLattice<T>::addCappedBondedIG()\n";
1295  CVarMPIBuffer param_buffer(m_comm);
1296 
1297  // get params
1298  param_buffer.receiveBroadcast(0);
1299  int tag=param_buffer.pop_int();
1300  string name = param_buffer.pop_string();
1301  double k=param_buffer.pop_double();
1302  double rbreak=param_buffer.pop_double();
1303  double maxforce=param_buffer.pop_double();
1304 
1305  console.XDebug()
1306  << "Got CappedBondedIG parameters: " << tag
1307  << " " << name.c_str() << " "
1308  << k << " " << rbreak << " " << maxforce << "\n";
1309  // setup InteractionGroup
1310  CCappedBondedIGP param;
1311  param.k=k;
1312  param.rbreak=rbreak;
1313  param.tag = tag;
1314  param.m_force_limit=maxforce;
1317 
1318  // set unbreakable if rbeak<0
1319  if(rbreak<0){
1320  B_PIS->setUnbreakable(true);
1321  console.XDebug() << "set bpis unbreakable\"n";
1322  }
1323  // recv broadcast connection data
1324  /*console.XDebug()
1325  << "rank=" << m_tml_comm.rank()
1326  << "TSubLattice<T>::addCappedBondedIG(): receiving conn_data.\n";
1327 
1328  vector<int> conn_data;
1329  m_tml_comm.recv_broadcast_cont(conn_data,0);
1330  console.XDebug()
1331  << "rank=" << m_tml_comm.rank()
1332  << "TSubLattice<T>::addBondedIG(): conn_data.size()="
1333  << conn_data.size() << "\n";
1334  */
1335  vector<int> vi(2,-1);
1336  for(size_t i=0;i<m_temp_conn[tag].size();i+=2)
1337  {
1338  vi[0] = (m_temp_conn[tag][i]);
1339  vi[1] = (m_temp_conn[tag][i+1]);
1340  B_PIS->tryInsert(vi);
1341  }
1342 
1343  // add InteractionGroup to map
1344  m_bpis.insert(make_pair(name,B_PIS));
1345 
1346  console.XDebug() << "end TSubLattice<T>::addCappedBondedIG()\n";
1347 }
1348 
1349 template <class T>
1351 {
1352  console.Error() << "TSubLattice<T>::addRotBondedIG() => trying to add rotational bonded IG to nonrotational model\n";
1353 }
1354 
1355 template <class T>
1357 {
1358  console.Error() << "TSubLattice<T>::addRotThermBondedIG() => trying to add rotational thermal bonded IG to nonrotational model\n";
1359 }
1360 
1365 template <class T>
1367 {
1368  console.XDebug() << "TSubLattice<T>::addShortBondedIG()\n";
1369  CVarMPIBuffer param_buffer(m_comm);
1370 
1371  // get params
1372  param_buffer.receiveBroadcast(0);
1373  int tag=param_buffer.pop_int();
1374  string name = param_buffer.pop_string();
1375  double k=param_buffer.pop_double();
1376  double rbreak=param_buffer.pop_double();
1377 
1378  console.XDebug()
1379  << "Got ShortBondedIG parameters: " << tag
1380  << " " << name.c_str() << " "
1381  << k << " " << rbreak << "\n";
1382  // setup InteractionGroup
1383  CBondedIGP param;
1384  param.k=k;
1385  param.rbreak=rbreak;
1386  param.tag = tag;
1389 
1390  // recv broadcast connection data
1391  /*console.XDebug()
1392  << "rank=" << m_tml_comm.rank()
1393  << "TSubLattice<T>::addShortBondedIG(): receiving conn_data.\n";
1394 
1395  vector<int> conn_data;
1396  m_tml_comm.recv_broadcast_cont(conn_data,0);
1397  console.XDebug()
1398  << "rank=" << m_tml_comm.rank()
1399  << "TSubLattice<T>::addShortBondedIG(): conn_data.size()="
1400  << conn_data.size() << "\n";*/
1401 
1402  vector<int> vi(2,-1);
1403  for(size_t i=0;i<m_temp_conn[param.tag].size();i+=2)
1404  {
1405  vi[0] = (m_temp_conn[param.tag][i]);
1406  vi[1] = (m_temp_conn[param.tag][i+1]);
1407  B_PIS->tryInsert(vi);
1408  }
1409 
1410  // add InteractionGroup to map
1411  m_bpis.insert(make_pair(name,B_PIS));
1412 
1413  console.XDebug() << "end TSubLattice<T>::addShortBondedIG()\n";
1414 }
1415 
1421 template <class T>
1423 {
1424  console.XDebug() << "TSubLattice<T>::addExIG()\n";
1425  CVarMPIBuffer pbuffer(m_comm);
1426 
1427  // get params
1428  pbuffer.receiveBroadcast(0);
1429  string s1 = pbuffer.pop_string();
1430  string s2 = pbuffer.pop_string();
1431 
1432  //console.XDebug()<< s1.c_str() << " " << s2.c_str() << "\n";
1433  map<string,AParallelInteractionStorage*>::iterator bonded_ig=m_bpis.find(s1);
1434  map<string,AParallelInteractionStorage*>::iterator dynamic_ig=m_dpis.find(s2);
1435  if((bonded_ig!=m_bpis.end())&&(dynamic_ig!=m_dpis.end()))
1436  {
1437  // map iterators point to [key,value] pairs, therefore it->second
1438  // is the pointer to the PIS here
1439  dynamic_ig->second->addExIG(bonded_ig->second);
1440  }
1441  else
1442  {
1443  console.Error() << "TSubLattice<T>::setExIG() - nonexisting interaction group \n";
1444  }
1445 
1446  console.XDebug() << "end TSubLattice<T>::addExIG()\n";
1447 }
1448 
1454 template <class T>
1456 {
1457  console.XDebug() << "TSubLattice<T>::removeIG()\n";
1458  CVarMPIBuffer pbuffer(m_comm);
1459 
1460  // get params
1461  pbuffer.receiveBroadcast(0);
1462  string igname = pbuffer.pop_string();
1463  map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.find(igname);
1464  if(iter!=m_dpis.end()){
1465  delete m_dpis[igname];
1466  m_dpis.erase(iter);
1467  } else {
1468  console.Error() << "TSubLattice<T>::removeIG() - nonexisting interaction group - ignore removal\n";
1469  }
1470 }
1471 
1472 
1476 template <class T>
1478 {
1479  console.XDebug() << "TSubLattice<T>::exchangePos() \n" ;
1480 
1481  m_ppa->exchange(&T::getExchangeValues,&T::setExchangeValues);
1482 
1483  console.XDebug() << "end TSubLattice<T>::exchangePos()\n" ;
1484 }
1485 
1489 template <class T>
1491 {
1492  console.XDebug() << "TSubLattice<T>::zeroForces()\n";
1493 
1494  // particles
1495  m_ppa->forAllParticles(&T::zeroForce);
1496  // trimeshes
1497  for(map<string,TriMesh*>::iterator iter=m_mesh.begin();
1498  iter!=m_mesh.end();
1499  iter++){
1500  (iter->second)->zeroForces();
1501  }
1502 
1503  // 2d meshes
1504  for(map<string,Mesh2D*>::iterator iter=m_mesh2d.begin();
1505  iter!=m_mesh2d.end();
1506  iter++){
1507  (iter->second)->zeroForces();
1508  }
1509  // walls
1510  for(typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.begin();iter!=m_WIG.end();iter++)
1511  {
1512  (iter->second)->zeroForce();
1513  }
1514  console.XDebug() << "end TSubLattice<T>::zeroForces() \n";
1515 }
1516 
1523 template <class T>
1525 {
1526  console.XDebug() << "TSubLattice<T>::calcForces() \n";
1527 
1528  // the walls
1529  for(typename map<string,AWallInteractionGroup<T>*>::iterator it=m_WIG.begin();it!=m_WIG.end();it++)
1530  {
1531  (it->second)->calcForces();
1532  }
1533  // single particle IGs
1534  for(
1535  typename NameIGroupMap::iterator siter=this->m_singleParticleInteractions.begin();
1536  siter != m_singleParticleInteractions.end();
1537  siter++
1538  )
1539  {
1540  (siter->second)->calcForces();
1541  }
1542  // dynamically created IGs
1543  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();iter!=m_dpis.end();iter++)
1544  {
1545  (iter->second)->calcForces();
1546  }
1547  // bonded IGs
1548  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();iter!=m_bpis.end();iter++)
1549  {
1550  (iter->second)->calcForces();
1551  }
1552  // last, but not least - damping
1553  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();iter!=m_damping.end();iter++)
1554  {
1555  (iter->second)->calcForces();
1556  }
1557 
1558  console.XDebug() << "end TSubLattice<T>::calcForces() \n";
1559 }
1560 
1567 template <class T>
1569 {
1570  m_dt = dt;
1571  console.XDebug() << "TSubLattice<T>::setTimeStepSize() \n";
1572 
1573  // the walls
1574  for(typename map<string,AWallInteractionGroup<T>*>::iterator it=m_WIG.begin();it!=m_WIG.end();it++)
1575  {
1576  (it->second)->setTimeStepSize(dt);
1577  }
1578  // single particle IGs
1579  for(
1580  typename NameIGroupMap::iterator siter=this->m_singleParticleInteractions.begin();
1581  siter != m_singleParticleInteractions.end();
1582  siter++
1583  )
1584  {
1585  (siter->second)->setTimeStepSize(dt);
1586  }
1587  // dynamically created IGs
1588  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();iter!=m_dpis.end();iter++)
1589  {
1590  (iter->second)->setTimeStepSize(dt);
1591  }
1592  // bonded IGs
1593  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();iter!=m_bpis.end();iter++)
1594  {
1595  (iter->second)->setTimeStepSize(dt);
1596  }
1597  // last, but not least - damping
1598  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();iter!=m_damping.end();iter++)
1599  {
1600  (iter->second)->setTimeStepSize(dt);
1601  }
1602 
1603  console.XDebug() << "end TSubLattice<T>::setTimeStepSize() \n";
1604 }
1605 
1611 template <class T>
1613 {
1614  console.XDebug() << "TSubLattice<T>::integrate \n";
1615  m_ppa->forAllParticles(&T::integrate,dt);
1616  m_ppa->forAllParticles(&T::rescale) ;
1617  console.XDebug() << "end TSubLattice<T>::integrate \n";
1618 }
1619 
1623 template <class T>
1625 {
1626  zeroForces();
1627  calcForces();
1628  integrate(m_dt);
1629 
1630  if (this->getParticleType() == "RotTherm")
1631  {
1632  this->oneStepTherm();
1633  }
1634 }
1635 
1639 template <class T>
1641 {
1642  zeroHeat(); // ???? combine?
1643  calcHeatFrict();
1644  calcHeatTrans();
1645  integrateTherm(m_dt);
1646  thermExpansion() ;
1647 }
1648 
1654 template <class T>
1656 {
1657  console.XDebug() << "TSubLattice<T>::integrateTherm \n";
1658  m_ppa->forAllParticles(&T::integrateTherm,dt);
1659 // m_ppa->forAllParticles(&T::rescale) ;
1660  console.XDebug() << "end TSubLattice<T>::integrateTherm \n";
1661 }
1662 
1663 template <class T>
1665 {
1666  console.XDebug() << "TSubLattice<T>::thermExpansion() \n";
1667  m_ppa->forAllParticles(&T::thermExpansion);
1668 // m_ppa->forAllParticles(&T::rescale) ;
1669  console.XDebug() << "end TSubLattice<T>::thermExpansion() \n";
1670 }
1671 
1675 template <class T>
1677 {
1678  console.XDebug() << "TSubLattice<T>::zeroHeat()\n";
1679 
1680  // particles
1681  m_ppa->forAllParticles(&T::zeroHeat);
1682 
1683 /*
1684 // walls
1685  for(typename vector<AWallInteractionGroup<T>*>::iterator iter=m_WIG.begin();iter!=m_WIG.end();iter++)
1686  {
1687  (*iter)->zeroForce();
1688  }
1689 */
1690  console.XDebug() << "end TSubLattice<T>::zeroHeat() \n";
1691 }
1692 
1699 template <class T>
1701 {
1702  console.XDebug() << "TSubLattice<T>::calcHeatFrict() \n";
1703 
1704  // dynamically created IGs
1705  for(
1706  typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1707  iter!=m_dpis.end();
1708  iter++
1709  )
1710  {
1711  (iter->second)->calcHeatFrict();
1712  }
1713 
1714  console.XDebug() << "end TSubLattice<T>::calcHeatFrict() \n";
1715 }
1716 
1717 template <class T>
1719 {
1720  console.XDebug() << "TSubLattice<T>::calcHeatTrans() \n";
1721 
1722 
1723  // dynamically created IGs
1724  for(
1725  typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1726  iter!=m_dpis.end();
1727  iter++
1728  )
1729  {
1730  (iter->second)->calcHeatTrans();
1731  }
1732  // bonded IGs
1733  for(
1734  typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();
1735  iter!=m_bpis.end();
1736  iter++
1737  )
1738  {
1739  (iter->second)->calcHeatTrans();
1740  }
1741 
1742  console.XDebug() << "end TSubLattice<T>::calcHeatTrans() \n";
1743 }
1744 
1748 template <class T>
1750 {
1751  m_ppa->rebuild();
1752 }
1753 
1757 template <class T>
1759 {
1760  CMPIBarrier barrier(m_worker_comm);
1761  m_pTimers->start("RebuildInteractions");
1762  m_pTimers->resume("NeighbourSearch");
1763  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();
1764  iter!=m_bpis.end();
1765  iter++)
1766  {
1767  console.Debug() << "exchg & rebuild BPIS " << iter->first << " at node " << m_rank << "\n";
1768  (iter->second)->exchange();
1769  (iter->second)->rebuild();
1770  }
1771 
1772  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1773  iter!=m_dpis.end();
1774  iter++)
1775  {
1776  console.Debug() << "exchg & rebuild DPIS " << iter->first << " at node " << m_rank << "\n";
1777  (iter->second)->exchange();
1778  m_pTimers->pause("RebuildInteractions");
1779  m_pTimers->pause("NeighbourSearch");
1780  barrier.wait("dpis::exchange");
1781  m_pTimers->resume("RebuildInteractions");
1782  m_pTimers->resume("NeighbourSearch");
1783  (iter->second)->rebuild();
1784  }
1785  resetDisplacements();
1786  m_pTimers->stop("RebuildInteractions");
1787 }
1788 
1792 template <class T>
1794 {
1795  console.Debug() << "CSubLattice<T>::searchNeighbors()\n";
1796  CMPIBarrier barrier(getWorkerComm());
1797  m_pTimers->start("NeighbourSearch");
1798  m_pTimers->start("RebuildParticleArray");
1799  rebuildParticleArray();
1800  m_pTimers->stop("RebuildParticleArray");
1801  m_pTimers->pause("NeighbourSearch");
1802  barrier.wait("PPA rebuild");
1803  rebuildInteractions();
1804  m_pTimers->stop("NeighbourSearch");
1805  console.Debug() << "end CSubLattice<T>::searchNeighbors()\n";
1806 }
1807 
1813 template <class T>
1815 {
1816  console.Debug() << "updateInteractions() \n";
1817  console.Debug() << "m_ppa->getTimeStamp() " << m_ppa->getTimeStamp() << " m_last_ns " << m_last_ns << "\n";
1818  bool need_update=false;
1819 
1820  m_pTimers->start("UpdateBondedInteractions");
1821  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();
1822  iter!=m_bpis.end();
1823  iter++)
1824  {
1825  bool n=(iter->second)->update();
1826  need_update=need_update || n;
1827  }
1828  m_pTimers->stop("UpdateBondedInteractions");
1829  if((m_ppa->getTimeStamp() > m_last_ns) || need_update)
1830  {
1831  m_pTimers->start("UpdateDynamicInteractions");
1832  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1833  iter!=m_dpis.end();
1834  iter++)
1835  {
1836  bool n=(iter->second)->update();
1837  need_update=need_update || n;
1838  }
1839  m_pTimers->stop("UpdateDynamicInteractions");
1840  for(typename map<string,AWallInteractionGroup<T>*>::iterator it=m_WIG.begin();
1841  it!=m_WIG.end();
1842  it++)
1843  {
1844  (it->second)->Update(m_ppa);
1845  }
1846  for(typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();
1847  iter!=m_damping.end();
1848  iter++){
1849  (iter->second)->update();
1850  }
1851  m_last_ns=m_ppa->getTimeStamp();
1852  }
1853 
1854  console.Debug() << "end TSubLattice<T>::updateInteractions()\n";
1855 }
1856 
1861 template <class T>
1863 {
1864  console.Debug() << "TSubLattice<T>::checkNeighbors()\n";
1865  CMPISGBufferLeaf buffer(m_comm,0,8);
1866  double mdsqr=0; // squared max. displacement
1867  double alpha=0.5*m_alpha;
1868  double srsqr=alpha*alpha; // squared search range
1869  vector<Vec3> displ; // displacements
1870  int res; // result 0/1
1871 
1872  // --- particles ---
1873  // get displacement data
1874  m_ppa->forAllParticlesGet(displ,&T::getDisplacement);
1875 
1876  // find maximum particle displacement
1877  vector<Vec3>::iterator it=displ.begin();
1878  while((it!=displ.end())&&(mdsqr<srsqr))
1879  {
1880  double sqdisp=(*it)*(*it);
1881  mdsqr = ((mdsqr < sqdisp) ? sqdisp : mdsqr);
1882  it++;
1883  //console.XDebug() << "sq. disp: " << sqdisp << "\n";
1884  }
1885  console.XDebug() << "max squared displacement " << mdsqr << "alpha^2 = " << srsqr << "\n";
1886  if (mdsqr>srsqr){
1887  res=1;
1888  } else {
1889  res=0;
1890  }
1891 
1892  // --- mesh ---
1893  // only needed if res==0
1894  if(res==0){
1895  for(map<string,TriMesh*>::iterator iter=m_mesh.begin();
1896  iter!=m_mesh.end();
1897  iter++){
1898  //std::cerr << "check mesh " << iter->first << std::endl;
1899  if(iter->second->hasMovedBy(alpha)){
1900  res=1;
1901  }
1902  }
1903  }
1904  buffer.append(res);
1905  buffer.send();
1906  console.Debug() << "end TSubLattice<T>::checkNeighbors()\n";
1907 }
1908 
1909 
1913 template <class T>
1915 {
1916  console.Debug() << "slave " << m_rank << " resetDisplacements()\n";
1917  m_ppa->forAllParticles(&T::resetDisplacement);
1918  for(map<string,TriMesh*>::iterator iter=m_mesh.begin();
1919  iter!=m_mesh.end();
1920  iter++){
1921  iter->second->resetCurrentDisplacement();
1922  }
1923  console.Debug() << "slave " << m_rank << " end resetDisplacements()\n";
1924 }
1925 
1930 template <class T>
1932 {
1933  console.Debug() << "TSubLattice<T>::moveParticleTo()\n";
1934  CVarMPIBuffer buffer(m_comm);
1935 
1936  buffer.receiveBroadcast(0); // get data from master
1937  int tag=buffer.pop_int();
1938  Vec3 mv=buffer.pop_vector();
1939  m_ppa->forParticleTag(tag,(void (T::*)(Vec3))(&T::moveToRel),mv);
1940  console.Debug() << "end TSubLattice<T>::moveParticleTo()\n";
1941 }
1942 
1947 template <class T>
1949 {
1950  console.Debug() << "TSubLattice<T>::moveTaggedParticlesBy()\n";
1951  CVarMPIBuffer buffer(m_comm);
1952 
1953  buffer.receiveBroadcast(0); // get data from master
1954  const int tag = buffer.pop_int();
1955  const Vec3 dx = buffer.pop_vector();
1956  m_ppa->forParticleTag(tag, (void (T::*)(Vec3))(&T::moveBy),dx);
1957  console.Debug() << "end TSubLattice<T>::moveTaggedParticlesBy()\n";
1958 }
1959 
1960 
1961 template <class T>
1962 void TSubLattice<T>::moveSingleParticleTo(int particleId, const Vec3 &posn)
1963 {
1964  m_ppa->forParticle(particleId, (void (T::*)(Vec3))(&T::moveTo), posn);
1965 }
1966 
1971 template <class T>
1973 {
1974  console.Debug() << "TSubLattice<T>::moveSingleNode()\n";
1975  CVarMPIBuffer buffer(m_comm);
1976 
1977  buffer.receiveBroadcast(0); // get data from master
1978  string name=string(buffer.pop_string());
1979  int id=buffer.pop_int();
1980  Vec3 disp=buffer.pop_vector();
1981 
1982  console.XDebug() << "name :" << name << " id : " << id << " disp " << disp << "\n";
1983 
1984  map<string,TriMesh*>::iterator tm=m_mesh.find(name);
1985  if (tm!=m_mesh.end()){
1986  (tm->second)->moveNode(id,disp);
1987  } else {
1988  map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(name);
1989  if(m2d!=m_mesh2d.end()){
1990  (m2d->second)->moveNode(id,disp);
1991  }
1992  }
1993  console.Debug() << "end TSubLattice<T>::moveSingleNode()\n";
1994 }
1995 
2000 template <class T>
2002 {
2003  console.Error() << "TSubLattice<T>::moveTaggedNodes() NOT IMPLEMENTED\n";
2004  throw
2005  std::runtime_error(
2006  "TSubLattice<T>::moveTaggedNodes() NOT IMPLEMENTED\n"
2007  );
2008 #if 0
2009  CVarMPIBuffer buffer(m_comm);
2010 
2011  buffer.receiveBroadcast(0); // get data from master
2012  string name=string(buffer.pop_string());
2013  int tag=buffer.pop_int();
2014  Vec3 disp=buffer.pop_vector();
2015 #endif
2016 }
2017 
2024 template <class T>
2025 void TSubLattice<T>::translateMeshBy(const std::string &meshName, const Vec3 &translation)
2026 {
2027  map<string,TriMesh*>::iterator tm=m_mesh.find(meshName);
2028  if (tm != m_mesh.end()){
2029  (tm->second)->translateBy(translation);
2030  }
2031  map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshName);
2032  if(m2d!=m_mesh2d.end()){
2033  (m2d->second)->translateBy(translation);
2034  }
2035 }
2036 
2037 template <class T>
2038 std::pair<double, int> TSubLattice<T>::findParticleNearestTo(const Vec3 &pt)
2039 {
2040  console.Debug() << "TSubLattice<T>::findParticleNearestTo: enter\n";
2041  const T *pClosest = NULL;
2042  double minDistSqrd = std::numeric_limits<double>::max();
2043 
2044  typename ParticleArray::ParticleIterator it =
2045  m_ppa->getInnerParticleIterator();
2046  while (it.hasNext())
2047  {
2048  const T &p = it.next();
2049  const double distSqrd = (pt - p.getPos()).norm2();
2050  if (distSqrd < minDistSqrd)
2051  {
2052  minDistSqrd = distSqrd;
2053  pClosest = &p;
2054  }
2055  }
2056  console.Debug() << "TSubLattice<T>::findParticleNearestTo: exit\n";
2057  return
2058  (
2059  (pClosest != NULL)
2060  ?
2061  std::make_pair(sqrt(minDistSqrd), pClosest->getID())
2062  :
2063  std::make_pair(std::numeric_limits<double>::max(), -1)
2064  );
2065 }
2066 
2070 template <class T>
2071 std::pair<int, Vec3> TSubLattice<T>::getParticlePosn(int particleId)
2072 {
2073  const T *particle = NULL;
2074  typename ParticleArray::ParticleIterator it =
2075  m_ppa->getInnerParticleIterator();
2076  while (it.hasNext())
2077  {
2078  const T &p = it.next();
2079  if (p.getID() == particleId)
2080  {
2081  particle = &p;
2082  }
2083  }
2084  if (particle != NULL)
2085  {
2086  return std::make_pair(particleId, particle->getPos());
2087  }
2088  return std::make_pair(-1,Vec3::ZERO);
2089 }
2090 
2094 template <class T>
2095 void TSubLattice<T>::getParticleData(const IdVector &particleIdVector)
2096 {
2097  console.Debug()
2098  << "TSubLattice<T>::getParticleData: enter\n";
2099  typedef std::set<int> IdSet;
2100  typedef std::vector<T> ParticleVector;
2101 
2102  ParticleVector particleVector;
2103  typename ParticleArray::ParticleIterator it =
2104  m_ppa->getInnerParticleIterator();
2105  if (particleIdVector.size() > 0)
2106  {
2107  IdSet idSet(particleIdVector.begin(), particleIdVector.end());
2108  console.Debug()
2109  << "TSubLattice<T>::getParticleData: iterating over particles\n";
2110  while (it.hasNext())
2111  {
2112  const T &p = it.next();
2113  if (idSet.find(p.getID()) != idSet.end())
2114  {
2115  particleVector.push_back(p);
2116  }
2117  }
2118  }
2119  else
2120  {
2121  m_ppa->getAllInnerParticles(particleVector);
2122  }
2123  console.Debug()
2124  << "TSubLattice<T>::getParticleData:"
2125  << " sending particle data of size " << particleVector.size() << "\n";
2126  m_tml_comm.send_gather_packed(particleVector, 0);
2127  console.Debug()
2128  << "TSubLattice<T>::getParticleData: exit\n";
2129 }
2130 
2134 template <class T>
2136 {
2137  console.Debug() << "TSubLattice<T>::tagParticleNearestTo()\n";
2138  CVarMPIBuffer buffer(m_comm);
2139 
2140  buffer.receiveBroadcast(0); // get data from master
2141  int tag=buffer.pop_int();
2142  int mask=buffer.pop_int();
2143  Vec3 pos=buffer.pop_vector();
2144 
2145  // warning - this is ugly
2146  T* part_ptr=m_ppa->getParticlePtrByPosition(pos);
2147  if(part_ptr!=NULL){
2148  int old_tag=part_ptr->getTag();
2149  int new_tag=(old_tag & (~mask)) | (tag & mask);
2150  part_ptr->setTag(new_tag);
2151 
2152  cout << "pos, realpos: " << pos << " " << part_ptr->getPos() << " old tag, new tag " << old_tag << " " << part_ptr->getTag() << endl;
2153  }
2154  console.Debug() << "end TSubLattice<T>::tagParticleNearestTo()\n";
2155 }
2156 
2161 template <class T>
2163 {
2164  console.Debug() << "TSubLattice<T>::setParticleNonDynamic()\n";
2165  CVarMPIBuffer buffer(m_comm);
2166 
2167  buffer.receiveBroadcast(0); // get data from master
2168  int tag=buffer.pop_int();
2169  m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamic));
2170  console.Debug() << "end TSubLattice<T>::setParticleNonDynamic()\n";
2171 }
2172 
2177 template <class T>
2179 {
2180  console.Debug() << "TSubLattice<T>::setParticleNonRot()\n";
2181  CVarMPIBuffer buffer(m_comm);
2182 
2183  buffer.receiveBroadcast(0); // get data from master
2184  int tag=buffer.pop_int();
2185  m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamicRot));
2186  console.Debug() << "end TSubLattice<T>::setParticleNonRot()\n";
2187 }
2188 
2193 template <class T>
2195 {
2196  console.Debug() << "TSubLattice<T>::setParticleNonTrans()\n";
2197  CVarMPIBuffer buffer(m_comm);
2198 
2199  buffer.receiveBroadcast(0); // get data from master
2200  int tag=buffer.pop_int();
2201  m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamicLinear));
2202  console.Debug() << "end TSubLattice<T>::setParticleNonTrans()\n";
2203 }
2204 
2208 template <class T>
2210 {
2211  console.Debug() << "TSubLattice<T>::setTaggedParticleVel()\n";
2212  CVarMPIBuffer buffer(m_comm);
2213 
2214  buffer.receiveBroadcast(0); // get data from master
2215  int tag=buffer.pop_int();
2216  Vec3 v=buffer.pop_vector();
2217  m_ppa->forParticleTag(tag,(void (T::*)(Vec3))(&T::setVel),v);
2218  console.XDebug() << "end TSubLattice<T>::setTaggedParticleVel()\n";
2219 }
2220 
2224 template <class T>
2226 {
2227  console.XDebug() << "TSubLattice<T>::moveWallBy()\n";
2228  CVarMPIBuffer buffer(m_comm);
2229 
2230  buffer.receiveBroadcast(0); // get data from master
2231  string wname=buffer.pop_string();
2232  Vec3 mv=buffer.pop_vector();
2233  typename map<string,CWall*>::iterator iter=m_walls.find(wname);
2234  if(iter!=m_walls.end())
2235  {
2236  (iter->second)->moveBy(mv);
2237  }
2238 }
2239 
2243 template <class T>
2245 {
2246  console.XDebug() << "TSubLattice<T>::setWallNormal()\n";
2247  CVarMPIBuffer buffer(m_comm);
2248 
2249  buffer.receiveBroadcast(0); // get data from master
2250  string wname=buffer.pop_string();
2251  Vec3 wn=buffer.pop_vector();
2252  typename map<string,CWall*>::iterator iter=m_walls.find(wname);
2253  if(iter!=m_walls.end())
2254  {
2255  (iter->second)->setNormal(wn);
2256  }
2257 }
2258 
2262 template <class T>
2264 {
2265  console.XDebug() << "TSubLattice<T>::applyForceToWall()\n";
2266  CVarMPIBuffer buffer(m_comm);
2267 
2268  buffer.receiveBroadcast(0); // get data from master
2269  string wname=buffer.pop_string();
2270  Vec3 f=buffer.pop_vector();
2271  typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.find(wname);
2272  if(iter!=m_WIG.end())
2273  {
2274  (iter->second)->applyForce(f);
2275  }
2276 }
2277 
2282 template <class T>
2284 {
2285  console.XDebug() << "TSubLattice<T>::setVelocityOfWall()\n";
2286  CVarMPIBuffer buffer(m_comm);
2287 
2288  buffer.receiveBroadcast(0); // get data from master
2289  string wname=buffer.pop_string();
2290  Vec3 v=buffer.pop_vector();
2291  typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.find(wname);
2292  if(iter!=m_WIG.end())
2293  {
2294  (iter->second)->setVelocity(v);
2295  }
2296 }
2297 
2301 template <class T>
2303 {
2304  console.Debug() << "TSubLattice<T>::setParticleVelocity()\n";
2305  CVarMPIBuffer buffer(m_comm);
2306 
2307  buffer.receiveBroadcast(0); // get data from master
2308  int id=buffer.pop_int();
2309  Vec3 mv=buffer.pop_vector();
2310  m_ppa->forParticle(id,(void (T::*)(Vec3))(&T::setVel),mv);
2311  console.XDebug() << "end TSubLattice<T>::setParticleVelocity()\n";
2312 }
2313 
2317 template <class T>
2319 {
2320  console.Debug() << "TSubLattice<T>::setParticleDensity()\n";
2321  CVarMPIBuffer buffer(m_comm);
2322 
2323  buffer.receiveBroadcast(0); // get data from master
2324  int tag=buffer.pop_int();
2325  int mask=buffer.pop_int();
2326  double rho=buffer.pop_double();
2327  m_ppa->forParticleTagMask(tag,mask,(void (T::*)(double))(&T::setDensity),rho);
2328  console.XDebug() << "end TSubLattice<T>::setParticleVelocity()\n";
2329 }
2330 
2331 
2335 template <class T>
2337 {
2338  console.Debug() << "TSubLattice<T>::sendDataToMaster()\n";
2339  vector<Vec3> positions;
2340  vector<double> radii;
2341  vector<int> ids;
2342 
2343  m_ppa->forAllParticlesGet(positions,(Vec3 (T::*)() const)(&T::getPos));
2344  m_ppa->forAllParticlesGet(radii,(double (T::*)() const)(&T::getRad));
2345  m_ppa->forAllParticlesGet(ids,(int (T::*)() const)(&T::getID));
2346 
2347  m_tml_comm.send_gather(positions,0);
2348  m_tml_comm.send_gather(radii,0);
2349  m_tml_comm.send_gather(ids,0);
2350 
2351  console.Debug() << "end TSubLattice<T>::sendDataToMaster()\n";
2352 }
2353 
2357 template <class T>
2359 {
2360  console.Debug()<<"TSubLattice<T>::countParticles()\n";
2361  CMPIVarSGBufferLeaf buffer(m_comm,0);
2362  //-- pack particles into buffer
2363  buffer.append(m_ppa->size());
2364  // send
2365  buffer.send();
2366 }
2367 
2371 template <class T>
2373 {
2374  cout<< "My Rank : " << m_rank << "\n" ;
2375  if(m_ppa!=NULL)
2376  {
2377  cout << *m_ppa << endl;
2378  }
2379 }
2380 
2381 template <class T>
2383 {
2384  cout << "Data: my rank : " << m_rank << "particles : \n" ;
2385  m_ppa->forAllParticles((void (T::*)())(&T::print));
2386 }
2387 
2388 template <class T>
2390 {
2391  console.Debug() << "time spent calculating force : " << forcetime << " sec\n";
2392  console.Debug() << "time spent communicating : " << commtime << " sec\n";
2393  console.Debug() << "time spent packing : " << packtime << " sec\n";
2394  console.Debug() << "time spent unpacking : " << unpacktime << " sec\n";
2395 }
2396 
2397 //-------------- FIELD FUNCTIONS ----------------
2398 
2399 
2400 
2404 template <class T>
2406 {
2407  // cout << "TSubLattice<T>::addScalarParticleField\n";
2408  string fieldname;
2409  int id,is_tagged;
2410 
2411  m_tml_comm.recv_broadcast_cont(fieldname,0);
2412  //cout << "recvd. fieldname: " << fieldname << "\n";
2413  m_tml_comm.recv_broadcast(id,0);
2414  //cout << "recvd. id: " << id << "\n";
2415  m_tml_comm.recv_broadcast(is_tagged,0);
2416  //cout << "recvd. is_tagged: " << is_tagged << "\n";
2417 
2418  typename T::ScalarFieldFunction rdf=T::getScalarFieldFunction(fieldname);
2419  ScalarParticleFieldSlave<T> *new_spfs;
2420  if(is_tagged==0)
2421  {
2422  new_spfs=new ScalarParticleFieldSlave<T>(&m_tml_comm,m_ppa,rdf);
2423  }
2424  else
2425  {
2426  int tag,mask;
2427  m_tml_comm.recv_broadcast(tag,0);
2428  console.XDebug() << "recvd. tag: " << tag << "\n";
2429  m_tml_comm.recv_broadcast(mask,0);
2430  console.XDebug() << "recvd. mask: " << mask << "\n";
2431  new_spfs=new ScalarParticleFieldSlaveTagged<T>(&m_tml_comm,m_ppa,rdf,tag,mask);
2432  }
2433  m_field_slaves.insert(make_pair(id,new_spfs));
2434 }
2435 
2439 template <class T>
2441 {
2442  console.XDebug() << "TSubLattice<T>::addVectorParticleField\n";
2443  string fieldname;
2444  int id,is_tagged;
2445 
2446  m_tml_comm.recv_broadcast_cont(fieldname,0);
2447  console.XDebug() << "recvd. fieldname: " << fieldname << "\n";
2448  m_tml_comm.recv_broadcast(id,0);
2449  console.XDebug() << "recvd. id: " << id << "\n";
2450  m_tml_comm.recv_broadcast(is_tagged,0);
2451  console.XDebug() << "recvd. is_tagged: " << is_tagged << "\n";
2452 
2453  typename T::VectorFieldFunction rdf=T::getVectorFieldFunction(fieldname);
2454  VectorParticleFieldSlave<T> *new_vpfs;
2455  if(is_tagged==0)
2456  {
2457  new_vpfs=new VectorParticleFieldSlave<T>(&m_tml_comm,m_ppa,rdf);
2458  }
2459  else
2460  {
2461  int tag,mask;
2462  m_tml_comm.recv_broadcast(tag,0);
2463  console.XDebug() << "recvd. tag: " << tag << "\n";
2464  m_tml_comm.recv_broadcast(mask,0);
2465  console.XDebug() << "recvd. mask: " << mask << "\n";
2466  new_vpfs=new VectorParticleFieldSlaveTagged<T>(&m_tml_comm,m_ppa,rdf,tag,mask);
2467  }
2468  m_field_slaves.insert(make_pair(id,new_vpfs));
2469 
2470  console.Debug() << "end TSubLattice<T>::addVectorParticleField\n";
2471 }
2472 
2473 
2477 template <class T>
2479 {
2480  console.XDebug() << "TSubLattice<T>::addScalarInteractionField\n";
2481  string fieldname;
2482  string igname;
2483  string igtype;
2484  int id,is_tagged,tag,mask,is_checked;
2485 
2486  m_tml_comm.recv_broadcast_cont(fieldname,0);
2487  console.XDebug() << "recvd. fieldname: " << fieldname << "\n";
2488  m_tml_comm.recv_broadcast(id,0);
2489  console.XDebug() << "recvd. id: " << id << "\n";
2490  m_tml_comm.recv_broadcast_cont(igname,0);
2491  console.XDebug() << "recvd. interaction group name: " << igname << "\n";
2492  m_tml_comm.recv_broadcast_cont(igtype,0);
2493  console.XDebug() << "recvd. interaction group name: " << igtype << "\n";
2494  m_tml_comm.recv_broadcast(is_tagged,0);
2495  console.XDebug() << "recvd. is_tagged: " << is_tagged << "\n";
2496 
2497  // get interaction group
2498  map<string,AParallelInteractionStorage*>::iterator it=m_dpis.find(igname);
2499  if(is_tagged==1)
2500  {
2501  m_tml_comm.recv_broadcast(tag,0);
2502  m_tml_comm.recv_broadcast(mask,0);
2503  }
2504  m_tml_comm.recv_broadcast(is_checked,0);
2505  console.XDebug() << "recvd. is_checked: " << is_checked << "\n";
2506 
2507  if(it!=m_dpis.end())
2508  {
2509  console.XDebug() << "found interaction group in dynamic\n";
2510  AFieldSlave* new_sifs;
2511  new_sifs=(it->second)->generateNewScalarFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2512  m_field_slaves.insert(make_pair(id,new_sifs));
2513  }
2514  else
2515  {
2516  it=m_bpis.find(igname);
2517  if(it!=m_bpis.end()){
2518  console.XDebug() << "found interaction group in bonded\n";
2519  AFieldSlave* new_sifs;
2520  new_sifs=(it->second)->generateNewScalarFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2521  m_field_slaves.insert(make_pair(id,new_sifs));
2522  }
2523  else // not in dynamic or bonded -> try damping
2524  {
2525  //typename map<string,CDampingGroup<T>*>::iterator it2=m_damping.find(igname);
2526  it=m_damping.find(igname);
2527  if(it!=m_damping.end()) // found it in damping
2528  {
2529  AFieldSlave* new_sifs;
2530  new_sifs=(it->second)->generateNewScalarFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2531  m_field_slaves.insert(make_pair(id,new_sifs));
2532  }
2533  else // still not found -> unknown name -> error
2534  {
2535  cerr << "ERROR : unknown interaction group name in addScalarInteractionField " << endl;
2536  }
2537  }
2538  }
2539 
2540  console.XDebug() << "end TSubLattice<T>::addScalarInteractionField\n";
2541 }
2542 
2546 template <class T>
2548 {
2549  console.Debug() << "TSubLattice<T>::addVectorInteractionField\n";
2550  string fieldname;
2551  string igname;
2552  string igtype;
2553  int id,is_tagged,tag,mask,is_checked;
2554 
2555  m_tml_comm.recv_broadcast_cont(fieldname,0);
2556  console.XDebug() << "recvd. fieldname: " << fieldname << "\n";
2557  m_tml_comm.recv_broadcast(id,0);
2558  console.XDebug() << "recvd. id: " << id << "\n";
2559  m_tml_comm.recv_broadcast_cont(igname,0);
2560  console.XDebug() << "recvd. interaction group name: " << igname << "\n";
2561  m_tml_comm.recv_broadcast_cont(igtype,0);
2562  console.XDebug() << "recvd. interaction group type: " << igtype << "\n";
2563  m_tml_comm.recv_broadcast(is_tagged,0);
2564  console.XDebug() << "recvd. is_tagged: " << is_tagged << "\n";
2565 
2566  // get interaction group
2567  map<string,AParallelInteractionStorage*>::iterator it=m_dpis.find(igname);
2568  if(is_tagged==1)
2569  {
2570  m_tml_comm.recv_broadcast(tag,0);
2571  m_tml_comm.recv_broadcast(mask,0);
2572  }
2573  m_tml_comm.recv_broadcast(is_checked,0);
2574  console.XDebug() << "recvd. is_checked: " << is_checked << "\n";
2575 
2576  if(it!=m_dpis.end())
2577  {
2578  console.XDebug() << "found interaction group in dynamic\n";
2579  AFieldSlave* new_sifs;
2580  new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2581  if(new_sifs!=NULL){
2582  m_field_slaves.insert(make_pair(id,new_sifs));
2583  } else {
2584  console.Error()<<"ERROR: could not generate Field Slave for field " << fieldname << "\n";
2585  }
2586  }
2587  else
2588  {
2589  it=m_bpis.find(igname);
2590  if(it!=m_bpis.end()){
2591  console.XDebug() << "found interaction group in bonded\n";
2592  AFieldSlave* new_sifs;
2593  new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2594  m_field_slaves.insert(make_pair(id,new_sifs));
2595  }
2596  else // not in dynamic or bonded -> try damping
2597  {
2598  //typename map<string,CDampingGroup<T>*>::iterator it2=m_damping.find(igname);
2599  it=m_damping.find(igname);
2600  if(it!=m_damping.end()) // found it in damping
2601  {
2602  AFieldSlave* new_sifs;
2603  new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2604  m_field_slaves.insert(make_pair(id,new_sifs));
2605  }
2606  else // still not found -> unknown name -> error
2607  {
2608  cerr << "ERROR : unknown interaction group name in addScalarInteractionField " << endl;
2609  }
2610  }
2611  }
2612 
2613  console.Debug() << "end TSubLattice<T>::addVectorInteractionField\n";
2614 }
2615 
2620 template <class T>
2622 {
2623  console.Debug() << "TSubLattice<T>::addVectorTriangleField()\n";
2624  string fieldname;
2625  string meshname;
2626  int id;
2627 
2628  // receive params
2629  m_tml_comm.recv_broadcast_cont(fieldname,0);
2630  console.XDebug() << "recvd. fieldname: " << fieldname << "\n";
2631  m_tml_comm.recv_broadcast_cont(meshname,0);
2632  console.XDebug() << "recvd. meshname: " << meshname << "\n";
2633  m_tml_comm.recv_broadcast(id,0);
2634  console.XDebug() << "recvd. id: " << id << "\n";
2635 
2636  map<string,TriMesh*>::iterator tm=m_mesh.find(meshname);
2637  // if meshname is in trimesh map
2638  if (tm!=m_mesh.end()){
2639  // get reader function
2640  Triangle::VectorFieldFunction rdf=Triangle::getVectorFieldFunction(fieldname);
2641  // check it
2642  if(rdf==NULL){
2643  console.Critical() << "NULL rdf !!!\n";
2644  }
2645  VectorTriangleFieldSlave* new_vfs=new VectorTriangleFieldSlave(&m_tml_comm,tm->second,rdf);
2646  m_field_slaves.insert(make_pair(id,new_vfs));
2647  } else {
2648  map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname);
2649  if(m2d!=m_mesh2d.end()){
2650  Edge2D::VectorFieldFunction rdf=Edge2D::getVectorFieldFunction(fieldname);
2651  // check it
2652  if(rdf==NULL){
2653  console.Critical() << "NULL rdf !!!\n";
2654  }
2655  VectorEdge2DFieldSlave* new_efs= new VectorEdge2DFieldSlave(&m_tml_comm,m2d->second,rdf);
2656  m_field_slaves.insert(make_pair(id,new_efs));
2657  }
2658  }
2659  console.Debug() << "end TSubLattice<T>::addVectorTriangleField()\n";
2660 }
2661 
2665 template <class T>
2667 {
2668  console.Debug() << "TSubLattice<T>::addScalarTriangleField()\n";
2669  string fieldname;
2670  string meshname;
2671  int id;
2672 
2673  // receive params
2674  m_tml_comm.recv_broadcast_cont(fieldname,0);
2675  console.XDebug() << "recvd. fieldname: " << fieldname << "\n";
2676  m_tml_comm.recv_broadcast_cont(meshname,0);
2677  console.XDebug() << "recvd. meshname: " << meshname << "\n";
2678  m_tml_comm.recv_broadcast(id,0);
2679  console.XDebug() << "recvd. id: " << id << "\n";
2680 
2681  // get reader function
2682  Triangle::ScalarFieldFunction rdf=Triangle::getScalarFieldFunction(fieldname);
2683  // check it
2684  if(rdf==NULL){
2685  console.Critical() << "NULL rdf !!!\n";
2686  }
2687  ScalarTriangleFieldSlave* new_vtfs=new ScalarTriangleFieldSlave(&m_tml_comm,m_mesh[meshname],rdf);
2688  m_field_slaves.insert(make_pair(id,new_vtfs));
2689  console.Debug() << "end TSubLattice<T>::addScalarTriangleField()\n";
2690 }
2691 
2695 template <class T>
2697 {
2698  console.XDebug() << "begin TSubLattice<T>::addVectorWallField()\n";
2699  string fieldname;
2700  string tmp_wallname;
2701  vector<string> wallnames;
2702  int nwall;
2703  int id;
2704 
2705  // receive params
2706  m_tml_comm.recv_broadcast_cont(fieldname,0);
2707  console.XDebug() << "recvd. fieldname: " << fieldname << "\n";
2708  m_tml_comm.recv_broadcast(nwall,0);
2709  console.XDebug() << "recvd. nwall: " << nwall << "\n";
2710  for(int i=0;i<nwall;i++){
2711  m_tml_comm.recv_broadcast_cont(tmp_wallname,0);
2712  wallnames.push_back(tmp_wallname);
2713  console.XDebug() << "recvd. wallname: " << tmp_wallname << "\n";
2714  tmp_wallname.clear();
2715  }
2716  m_tml_comm.recv_broadcast(id,0);
2717  console.XDebug() << "recvd. id: " << id << "\n";
2718 
2719  // check validity of 1st wall name
2720  map<string,CWall*>::iterator cwalliter=m_walls.find(*(wallnames.begin()));
2721  if(cwalliter==m_walls.end()){ // 1st wall name invalid -> exit
2722  std::stringstream msg;
2723  msg
2724  << "ERROR in addVectorWallField: wallname '"
2725  << *(wallnames.begin()) << " 'invalid. Valid wall names: ";
2726  for (map<string,CWall*>::const_iterator it = m_walls.begin(); it != m_walls.end(); it++)
2727  {
2728  msg << "'" << it->first << "' ";
2729  }
2730  console.Error() << msg.str() << "\n";
2731  throw std::runtime_error(msg.str());
2732  } else { // first wall valid -> try to get slave
2733  // get summation flag from wall
2734  int sumflag=(cwalliter->second)->getFieldSummationFlag(fieldname);
2735  // if process 1, send summation flag back to master
2736  if(m_tml_comm.rank()==1){
2737  m_tml_comm.send(sumflag,0);
2738  }
2739  m_tml_comm.barrier();
2740 
2741  //field slave
2742  AWallFieldSlave* new_fs=(cwalliter->second)->generateVectorFieldSlave(&m_tml_comm,fieldname);
2743 
2744  // try to insert other walls
2745  vector<string>::iterator niter=wallnames.begin();
2746  if(niter!=wallnames.end()) niter++ ; // jump past 1st wall - already got it
2747  while(niter!=wallnames.end()){
2748  string wname=*niter;
2749  map<string,CWall*>::iterator iter=m_walls.find(wname);
2750  if(iter==m_walls.end()){ // wall name invalid -> exit
2751  std::stringstream msg;
2752  msg
2753  << "ERROR in addVectorWallField: wallname '"
2754  << wname << " 'invalid";
2755  for (map<string,CWall*>::const_iterator it = m_walls.begin(); it != m_walls.end(); it++)
2756  {
2757  msg << "'" << it->first << "' ";
2758  }
2759 
2760  console.Error() << msg.str() << "\n";
2761  throw std::runtime_error(msg.str());
2762  } else {
2763  new_fs->addWall(iter->second);
2764  }
2765  niter++;
2766  }
2767  if(new_fs!=NULL){
2768  m_field_slaves.insert(make_pair(id,new_fs));
2769  } else {
2770  console.Error() << "ERROR in addVectorWallField: got NULL Slave\n";
2771  }
2772  }
2773 
2774  console.XDebug() << "end TSubLattice<T>::addVectorWallField()\n";
2775 }
2776 
2780 template <class T>
2782 {
2783  console.Debug() << "TSubLattice<T>::sendFieldData()\n";
2784  // receive id's of field to send
2785  int id;
2786  m_tml_comm.recv_broadcast(id,0);
2787  console.Debug() << "received field id " << id << " for data collection\n" ;
2788  if(m_field_slaves[id] != NULL)
2789  {
2790  m_field_slaves[id]->sendData();
2791  }
2792  else
2793  { // can not happen
2794  cerr << "NULL pointer in m_field_slaves!" << endl;
2795  }
2796  // call the sending function of the field
2797  console.Debug() << "end TSubLattice<T>::sendFieldData()\n";
2798 }
2799 
2800 
2801 // ---- Checkpointing ----------
2805 template <class T>
2806 void TSubLattice<T>::saveSnapShotData(std::ostream &oStream)
2807 {
2808  // get precision of output stream and set it to 9 significant digits
2809  std::streamsize oldprec=oStream.precision(9);
2810 
2811  //
2812  // Save particle check-point data
2813  //
2814  ParticleArray &particleArray = dynamic_cast<ParticleArray &>(*m_ppa);
2816  particleIt(particleArray.getInnerParticleIterator());
2817 
2818  const std::string delim = "\n";
2819 
2820  oStream << particleIt.getNumRemaining() << delim;
2821  while (particleIt.hasNext()) {
2822  particleIt.next().saveSnapShotData(oStream);
2823  oStream << delim;
2824  }
2825 
2826  //
2827  // Save Bonded interaction check-point data.
2828  //
2829  typedef std::map<string,AParallelInteractionStorage*> NameBondedInteractionsMap;
2830  typename NameBondedInteractionsMap::iterator it;
2831  oStream << m_bpis.size() << delim;
2832  for (it = m_bpis.begin(); it != m_bpis.end(); it++) {
2833  it->second->saveSnapShotData(oStream);
2834  oStream << delim;
2835  }
2836 
2837  // dump trimeshdata (if exists)
2838  oStream << "TMIG " << m_mesh.size() << delim;
2839  for(typename map<string,TriMesh*>::iterator tm_iter=m_mesh.begin();
2840  tm_iter!=m_mesh.end();
2841  tm_iter++){
2842  oStream << tm_iter->first << delim;
2843  tm_iter->second->writeCheckPoint(oStream,delim);
2844  }
2845 
2846  // restore output stream to old precision
2847  oStream.precision(oldprec);
2848 }
2849 
2853 template <class T>
2854 void TSubLattice<T>::saveCheckPointData(std::ostream &oStream)
2855 {
2856  // get precision of output stream and set it to 12 significant digits
2857  std::streamsize oldprec=oStream.precision(16);
2858 
2859  const std::string delim = "\n";
2860 
2861  //
2862  // Save particle check-point data
2863  //
2864  m_ppa->saveCheckPointData(oStream);
2865 
2866  //
2867  // Save Bonded interaction check-point data.
2868  //
2869  typedef std::map<string,AParallelInteractionStorage*> NameBondedInteractionsMap;
2870  typename NameBondedInteractionsMap::iterator it;
2871  oStream << m_bpis.size() << delim;
2872  for (it = m_bpis.begin(); it != m_bpis.end(); it++) {
2873  it->second->saveCheckPointData(oStream);
2874  oStream << delim;
2875  }
2876 
2877  //
2878  // Save Non-bonded interaction check-point data
2879  //
2880  int count_save=0;
2881  for(std::map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
2882  iter!=m_dpis.end();
2883  iter++){
2884  if(iter->second->willSave()) count_save++;
2885  }
2886  oStream << count_save << delim;
2887  for(std::map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
2888  iter!=m_dpis.end();
2889  iter++){
2890  if(iter->second->willSave()) {
2891  iter->second->saveCheckPointData(oStream);
2892  oStream << delim;
2893  }
2894  }
2895 
2896  // Save walls (name, pos, normal)
2897  oStream << "Walls " << m_walls.size() << delim;
2898  for(map<string,CWall*>::iterator w_iter=m_walls.begin();
2899  w_iter!=m_walls.end();
2900  w_iter++){
2901  oStream << w_iter->first << delim;
2902  w_iter->second->writeCheckPoint(oStream,delim);
2903  }
2904 
2905  // dump trimeshdata (if exists)
2906  oStream << "TriMesh " << m_mesh.size() << delim;
2907  for(typename map<string,TriMesh*>::iterator tm_iter=m_mesh.begin();
2908  tm_iter!=m_mesh.end();
2909  tm_iter++){
2910  oStream << tm_iter->first << delim;
2911  tm_iter->second->writeCheckPoint(oStream,delim);
2912  }
2913  // dump 2D mesh data (if exists)
2914  oStream << "Mesh2D " << m_mesh2d.size() << delim;
2915  for(typename map<string,Mesh2D*>::iterator tm_iter=m_mesh2d.begin();
2916  tm_iter!=m_mesh2d.end();
2917  tm_iter++){
2918  oStream << tm_iter->first << delim;
2919  tm_iter->second->writeCheckPoint(oStream,delim);
2920  }
2921  // restore output stream to old precision
2922  oStream.precision(oldprec);
2923 }
2924 
2925 template <class T>
2926 void TSubLattice<T>::loadCheckPointData(std::istream &iStream)
2927 {
2928  // get particles
2929  m_ppa->loadCheckPointData(iStream);
2930 
2931  // rebuild neighbor table
2932  CMPIBarrier barrier(getWorkerComm());
2933  m_ppa->rebuild();
2934  barrier.wait("PPA rebuild");
2935 
2936  //-- get bonds --
2937  // get nr. of bonded interaction group in the checkpoint file
2938  int nr_bonded_ig;
2939  iStream >> nr_bonded_ig;
2940 
2941  // compare with existing bonded particle interaction storage (bpis)
2942  // barf if not equal
2943  if(nr_bonded_ig!=m_bpis.size()){
2944  std::cerr << "number of bonded interaction groups differ between snapshot and script!" << std::endl;
2945  } else { // numbers fit -> read data
2946  for (map<string,AParallelInteractionStorage*>::iterator it = m_bpis.begin();
2947  it != m_bpis.end();
2948  it++) { // for all interaction groups
2949  it->second->loadCheckPointData(iStream);
2950  }
2951  }
2952  //-- get nonbonded interactions --
2953  // get nr. of non-bonded interaction group in the checkpoint file
2954  int nr_nonbonded_ig;
2955  iStream >> nr_nonbonded_ig;
2956 
2957  // compare with existing non-bonded (dynamic) particle interaction storage (dpis)
2958  // barf if not equal
2959  if(nr_nonbonded_ig!=m_dpis.size()){
2960  std::cerr << "number of dynamic interaction groups differ between snapshot and script!" << std::endl;
2961  } else { // numbers fit -> read data
2962  for (map<string,AParallelInteractionStorage*>::iterator it = m_dpis.begin();
2963  it != m_dpis.end();
2964  it++) { // for all interaction groups
2965  it->second->loadCheckPointData(iStream);
2966  }
2967  }
2968  //--- load walls ---
2969  string token;
2970  iStream >> token;
2971  if(token!="Walls") { // found wrong token -> barf
2972  std::cerr << "expected Walls , got " << token << std::endl;
2973  }
2974  // nr. of walls
2975  int nwalls;
2976  iStream >> nwalls;
2977  // read wall names & data
2978  string wname;
2979  for(int i=0;i<nwalls;i++){
2980  CWall* new_wall = new CWall();
2981  iStream >> wname;
2982  new_wall->loadCheckPoint(iStream);
2983  m_walls[wname]=new_wall;
2984  }
2985  // --- load meshes --
2986  int nmesh;
2987  string mname;
2988  // Trimesh (3D)
2989  iStream >> token;
2990  if(token!="TriMesh") { // found wrong token -> barf
2991  std::cerr << "expected TriMesh , got " << token << std::endl;
2992  }
2993  // nr. of meshes
2994  iStream >> nmesh;
2995  // read wall names & data
2996  for(int i=0;i<nmesh;i++){
2997  TriMesh* new_tm=new TriMesh();
2998  iStream >> mname;
2999  new_tm->loadCheckPoint(iStream);
3000  m_mesh.insert(make_pair(mname,new_tm));
3001  }
3002  // Mesh2D
3003  iStream >> token;
3004  if(token!="Mesh2D") { // found wrong token -> barf
3005  std::cerr << "expected Mesh2D , got " << token << std::endl;
3006  }
3007  // nr. of meshes
3008  iStream >> nmesh;
3009  // read wall names & data
3010  for(int i=0;i<nmesh;i++){
3011  Mesh2D* new_m2d=new Mesh2D();
3012  iStream >> mname;
3013  new_m2d->loadCheckPoint(iStream);
3014  m_mesh2d.insert(make_pair(mname,new_m2d));
3015  }
3016 }
3017 
3018 // -- mesh data exchange functions --
3019 
3023 template <class T>
3025 {
3026  console.XDebug() << "TSubLattice<T>::getMeshNodeRef()\n";
3027  vector<int> ref_vec;
3028 
3029  // MPI buffer
3030  CVarMPIBuffer param_buffer(m_comm);
3031  // receive mesh name
3032  param_buffer.receiveBroadcast(0);
3033  string meshname=param_buffer.pop_string();
3034  console.XDebug() << "Mesh name: " << meshname << "\n";
3035 
3036  // find mesh & collect node ids into array
3037  map<string,TriMesh*>::iterator tm=m_mesh.find(meshname);
3038  if (tm!=m_mesh.end()){
3039  for(TriMesh::corner_iterator iter=(tm->second)->corners_begin();
3040  iter!=(tm->second)->corners_end();
3041  iter++){
3042  ref_vec.push_back(iter->getID());
3043  }
3044  } else {
3045  map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname);
3046  if(m2d!=m_mesh2d.end()){
3047  for(Mesh2D::corner_iterator iter=(m2d->second)->corners_begin();
3048  iter!=(m2d->second)->corners_end();
3049  iter++){
3050  ref_vec.push_back(iter->getID());
3051  }
3052  } else {
3053  console.Critical() << "ERROR - WRONG MESH NAME IN getMeshNodeRef() !! \n";
3054  }
3055  }
3056  // send back to master
3057  m_tml_comm.send_gather(ref_vec,0);
3058 
3059  console.XDebug() << "end TSubLattice<T>::getMeshNodeRef()\n";
3060 }
3061 
3065 template <class T>
3067 {
3068  console.XDebug() << "TSubLattice<T>::getMeshFaceRef()\n";
3069  vector<int> ref_vec;
3070 
3071  // MPI buffer
3072  CVarMPIBuffer param_buffer(m_comm);
3073  // receive mesh name
3074  param_buffer.receiveBroadcast(0);
3075  string meshname=param_buffer.pop_string();
3076  console.XDebug() << "Mesh name: " << meshname << "\n";
3077 
3078  // find mesh & collect node ids into array
3079  map<string,TriMesh*>::iterator tm=m_mesh.find(meshname);
3080  if (tm!=m_mesh.end()){
3081  for(TriMesh::triangle_iterator iter=(tm->second)->triangles_begin();
3082  iter!=(tm->second)->triangles_end();
3083  iter++){
3084  ref_vec.push_back(iter->getID());
3085  }
3086  } else {
3087  map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname);
3088  if(m2d!=m_mesh2d.end()){
3089  for(Mesh2D::edge_iterator iter=(m2d->second)->edges_begin();
3090  iter!=(m2d->second)->edges_end();
3091  iter++){
3092  ref_vec.push_back(iter->getID());
3093  }
3094  } else {
3095  console.Critical() << "ERROR - WRONG MESH NAME IN getMeshNodeRef() !! \n";
3096  }
3097  }
3098  // send back to master
3099  m_tml_comm.send_gather(ref_vec,0);
3100 
3101  console.XDebug() << "end TSubLattice<T>::getMeshNodeRef()\n";
3102 }
3103 
3107 template <class T>
3109 {
3110  console.XDebug() << "TSubLattice<T>::getMesh2DStress()\n";
3111  // receive mesh name
3112  // MPI buffer
3113  CVarMPIBuffer param_buffer(m_comm);
3114  param_buffer.receiveBroadcast(0);
3115  string meshname=param_buffer.pop_string();
3116  console.XDebug() << "Mesh name: " << meshname << "\n";
3117 
3118  // find mesh & collect data
3119  map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname);
3120  if(m2d!=m_mesh2d.end()){
3121  vector<pair<int,Vec3> > data_vec;
3122  // get data
3123  data_vec=(m2d->second)->forAllEdgesGetIndexed(&Edge2D::getForceDensity);
3124 
3125  // send data to master
3126  m_tml_comm.send_gather(data_vec,0);
3127  } else {
3128  console.Critical() << "ERROR - WRONG MESH NAME IN getMesh2DStress() !! \n";
3129  }
3130 
3131  console.XDebug() << "end TSubLattice<T>::getMesh2DStress()\n";
3132 }
3133 
3137 template <class T>
3139 {
3140  console.XDebug() << "TSubLattice<T>::getTriMeshStress(): enter\n";
3141  // receive mesh name
3142  // MPI buffers
3143  CVarMPIBuffer param_buffer(m_comm);
3144  param_buffer.receiveBroadcast(0);
3145  const std::string meshName = param_buffer.pop_string();
3146  console.XDebug() << "Mesh name: " << meshName << "\n";
3147 
3148  // find mesh & collect data
3149  map<string,TriMesh*>::iterator m=m_mesh.find(meshName);
3150  if(m != m_mesh.end()){
3151  vector<pair<int,Vec3> > data_vec;
3152  // get data
3153  data_vec=(m->second)->forAllTrianglesGetIndexed(&Triangle::getForce);
3154 
3155  // send data to master
3156  m_tml_comm.send_gather(data_vec,0);
3157  } else {
3158  std::stringstream msg;
3159  msg << "Invalid mesh name: " << meshName << ". No such triangular mesh.";
3160  throw std::runtime_error(msg.str().c_str());
3161  }
3162 
3163  console.XDebug() << "TSubLattice<T>::getTriMeshStress(): exit\n";
3164 }