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"
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"
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"
78 #include "Model/BodyForceGroup.h"
87 #include <boost/limits.hpp>
88 using std::runtime_error;
118 m_singleParticleInteractions(),
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),
149 m_dims = param.processDims();
157 m_alpha=param.alpha();
158 m_nrange=param.nrange();
172 console.
Debug() <<
"TSubLattice<T>::~TSubLattice(): enter\n";
174 <<
"TSubLattice<T>::~TSubLattice():"
175 <<
" deleting wall interaction groups...\n";
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";
200 console.
XDebug() <<
"TSubLattice<T>::initNeighborTable(" << min <<
"," << max <<
")\n";
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);
209 Vec3 nmin=min-0.5*grow;
210 Vec3 nmax=max+0.5*grow;
211 console.
XDebug() <<
"range=" << m_nrange <<
", new min,max: " << nmin <<
", " << nmax <<
"\n";
220 console.
XDebug() <<
"end TSubLattice<T>::initNeighborTable\n";
226 T::setDo2dCalculations(do2d);
232 return m_ppa->getInnerSize();
245 console.
XDebug() <<
"TSubLattice<T>::initNeighborTable(" << min <<
"," << max <<
") circ\n";
246 double xsize,ysize,zsize;
251 xsize=max.X()-min.X();
252 if(fabs(xsize/m_nrange-lrint(xsize/m_nrange))>1e-6)
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";
263 xsize=max.X()-min.X();
264 xsize=m_nrange*ceil(xsize/m_nrange);
269 ysize=max.Y()-min.Y();
270 if(fabs(ysize/m_nrange-lrint(ysize/m_nrange))>1e-6)
272 console.
Critical() <<
"circular y-dimension doesn't fit range !\n";
278 ysize=max.Y()-min.Y();
279 ysize=m_nrange*ceil(ysize/m_nrange);
284 zsize=max.Z()-min.Z();
285 if(fabs(zsize/m_nrange-lrint(zsize/m_nrange))>1e-6)
287 console.
Critical() <<
"circular z-dimension doesn't fit range !\n";
293 zsize=max.Z()-min.Z();
294 zsize=m_nrange*ceil(zsize/m_nrange);
296 Vec3 grow=
Vec3(xsize,ysize,zsize)-(max-min);
297 Vec3 nmin=min-0.5*grow;
298 Vec3 nmax=max+0.5*grow;
299 console.
XDebug() <<
"range, new min, max: " << m_nrange <<
" " << nmin << nmax <<
"\n";
306 console.
XDebug() <<
"end TSubLattice<T>::initNeighborTable (circ)\n";
318 console.
XDebug() <<
"TSubLattice<T>::receiveParticles: enter\n";
320 vector<T> recv_buffer;
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);
327 barrier.
wait(
"TSubLattice<T>::receiveParticles");
329 console.
XDebug() <<
"TSubLattice<T>::receiveParticles: exit\n";
341 console.
XDebug() <<
"TSubLattice<T>::receiveConnections: enter\n";
343 vector<int> recv_buffer;
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)
351 if ( (m_ppa->getParticlePtrByIndex( *(it+1)) == NULL ) ||
352 (m_ppa->getParticlePtrByIndex( *(it+2)) == NULL ) )
356 m_temp_conn[*(it)].push_back(*(it+1));
357 m_temp_conn[*(it)].push_back(*(it+2));
360 barrier.
wait(
"TSubLattice<T>::receiveConnections");
362 console.
XDebug() <<
"TSubLattice<T>::receiveConnections: exit\n";
372 console.
XDebug() <<
"TSubLattice<T>::addWall: enter\n" ;
377 Vec3 ipos=param_buffer.pop_vector();
378 Vec3 inorm=param_buffer.pop_vector();
380 m_walls[name]=
new CWall(ipos,inorm);
382 console.
XDebug() <<
"TSubLattice<T>::addWall: exit\n" ;
391 console.
XDebug() <<
"TSubLattice<T>::addElasticWIG: enter\n" ;
397 CEWallIGP* wigp=extractEWallIGP(¶m_buffer);
399 string wallname=wigp->getWallName();
400 map<string,CWall*>::iterator iter=m_walls.find(wallname);
401 if(iter!=m_walls.end()){
408 newCEWIG->Update(m_ppa);
409 m_WIG.insert(make_pair(wigp->getName(),newCEWIG));
411 std::stringstream msg;
412 msg <<
"wall name '" << wallname <<
"' not found in map of walls";
413 throw std::runtime_error(msg.str().c_str());
417 console.
XDebug() <<
"TSubLattice<T>::addElasticWIG: exit\n" ;
427 console.
XDebug() <<
"TSubLattice<T>::addBondedWIG: enter\n" ;
433 CBWallIGP* wigp=extractBWallIGP(¶m_buffer);
435 string wallname=wigp->getWallName();
436 map<string,CWall*>::iterator iter=m_walls.find(wallname);
437 if(iter!=m_walls.end()){
439 newCBWIG->Update(m_ppa);
440 m_WIG.insert(make_pair(wigp->getName(),newCBWIG));
442 console.
Error() <<
"wall name " << wallname <<
" not found in map of walls\n";
446 console.
XDebug() <<
"TSubLattice<T>::addBondedWIG: exit\n" ;
455 console.
XDebug() <<
"TSubLattice<T>::addDirBondedWIG: enter\n" ;
463 string wallname=wigp->getWallName();
464 map<string,CWall*>::iterator iter=m_walls.find(wallname);
465 if(iter!=m_walls.end()){
467 newCDWIG->Update(m_ppa);
468 m_WIG.insert(make_pair(wigp->getName(),newCDWIG));
470 console.
Error() <<
"wall name " << wallname <<
" not found in map of walls\n";
474 console.
XDebug() <<
"TSubLattice<T>::addDirBondedWIG: exit\n" ;
483 console.
XDebug() <<
"TSubLattice<T>::getWallPosition: enter\n" ;
491 console.
XDebug() <<
"Wall name: " << wname <<
"\n";
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";
499 pos=
Vec3(0.0,0.0,0.0);
504 m_tml_comm.send_gather(vpos,0);
505 console.
XDebug() <<
"TSubLattice<T>::getWallPosition: exit\n" ;
514 console.
XDebug() <<
"TSubLattice<T>::getWallForce: enter\n" ;
522 console.
XDebug() <<
"Wall name: " << wname <<
"\n";
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";
530 force=
Vec3(0.0,0.0,0.0);
534 vforce.push_back(force);
535 m_tml_comm.send_gather(vforce,0);
536 console.
XDebug() <<
"TSubLattice<T>::getWallForce: exit\n" ;
545 console.
XDebug() <<
"TSubLattice<T>::addViscWIG: enter\n" ;
551 CVWallIGP* wigp=extractVWallIGP(¶m_buffer);
553 string wallname=wigp->getWallName();
554 map<string,CWall*>::iterator iter=m_walls.find(wallname);
555 if(iter!=m_walls.end()){
557 newCVWIG->Update(m_ppa);
558 m_WIG.insert(make_pair(wigp->getName(),newCVWIG));
560 console.
Error() <<
"wall name " << wallname <<
" not found in map of walls\n";
564 console.
XDebug() <<
"TSubLattice<T>::addViscWIG: exit\n" ;
573 console.
XDebug() <<
"TSubLattice<T>::addPairIG()\n";
579 console.
XDebug()<<
"PIG type: " << type.c_str() <<
"\n";
581 console.
XDebug()<<
"PIG name: " << name.c_str() <<
"\n";
583 doAddPIG(name,type,param_buffer,
false);
585 console.
XDebug() <<
"end TSubLattice<T>::addPairIG()\n";
594 console.
XDebug() <<
"TSubLattice<T>::addTaggedPairIG()\n";
600 console.
XDebug()<<
"PIG type: " << type.c_str() <<
"\n";
602 console.
XDebug()<<
"PIG name: " << name.c_str() <<
"\n";
604 doAddPIG(name,type,param_buffer,
true);
606 console.
XDebug() <<
"end TSubLattice<T>::addTaggedPairIG()\n";
622 if(type==
"Elastic") {
625 eigp.m_scaling=
static_cast<bool>(param_buffer.
pop_int());
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";
639 }
else if (type==
"Friction") {
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";
650 }
else if (type==
"AdhesiveFriction") {
658 <<
"k,mu,k_s,dt,r_cut: " << figp.k <<
" , " << figp.mu <<
" , "
659 << figp.k_s <<
" , " << figp.dt <<
" " << figp.r_cut <<
"\n";
662 }
else if (type==
"FractalFriction") {
668 console.
XDebug() <<
"k,mu_0,k_s,dt: " << figp.k <<
" , " << figp.mu_0 <<
" , "
669 << figp.k_s <<
" , " << figp.dt <<
"\n";
674 figp.nx=param_buffer.
pop_int();
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]);
683 for(
int i=0;i<figp.nx*figp.
ny;i++)
690 }
else if(type==
"VWFriction") {
699 <<
"k,mu,k_s,dt,alpha: " << figp.k <<
" , " << figp.mu <<
" , "
700 << figp.k_s <<
" , " << figp.dt <<
"\n";
703 }
else if(type==
"RotElastic"){
708 }
else if (type==
"RotFriction"){
715 rfigp.scaling=
static_cast<bool>(param_buffer.
pop_int());
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";
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";
733 }
else if (type ==
"RotThermElastic") {
738 <<
"k=" << eigp.m_kr <<
" , "
739 <<
"diffusivity=" << eigp.diffusivity <<
"\n";
750 }
else if (type ==
"RotThermFriction") {
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";
775 }
else if(type==
"HertzianElastic") {
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();
789 }
else if(type==
"HertzianViscoElasticFriction") {
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();
807 }
else if(type==
"HertzianViscoElastic") {
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();
822 }
else if(type==
"LinearDashpot") {
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();
837 cerr <<
"Unknown interaction group name "
839 <<
" in TSubLattice<T>::addPairIG()" << endl;
843 if(res) m_dpis.insert(make_pair(name,new_pis));
855 console.
XDebug() <<
"TSubLattice<T>::addTriMesh()\n";
860 vector<MeshNodeData> node_recv_buffer;
861 vector<MeshTriData> tri_recv_buffer;
868 console.
XDebug()<<
"TriMesh name: " << name.c_str() <<
"\n";
871 m_tml_comm.recv_broadcast_cont_packed(node_recv_buffer,0);
872 console.
XDebug() <<
"recvd " << node_recv_buffer.size() <<
" nodes \n";
875 m_tml_comm.recv_broadcast_cont_packed(tri_recv_buffer,0);
876 console.
XDebug() <<
"recvd " << tri_recv_buffer.size() <<
" triangles \n";
880 new_tm->
LoadMesh(node_recv_buffer,tri_recv_buffer);
882 m_mesh.insert(make_pair(name,new_tm));
884 console.
XDebug() <<
"end TSubLattice<T>::addTriMesh()\n";
893 console.
XDebug() <<
"TSubLattice<T>::addTriMeshIG()\n";
903 console.
XDebug()<<
"TriMeshIG type: " << type.c_str() <<
"\n";
905 console.
XDebug()<<
"TriMeshIG name: " << name.c_str() <<
"\n";
907 console.
XDebug()<<
"TriMeshIG mesh name: " << meshname.c_str() <<
"\n";
911 if (m_mesh.find(meshname) != m_mesh.end())
913 tmp = m_mesh[meshname];
916 throw runtime_error(
"unknown mesh name in TSubLattice<T>::addTriMeshIG:" + meshname);
925 m_dpis.insert(make_pair(name,new_pis));
927 throw runtime_error(
"unknown type in TSubLattice<T>::addTriMeshIG:" + type);
931 console.
XDebug() <<
"end TSubLattice<T>::addTriMeshIG()\n";
940 console.
XDebug() <<
"TSubLattice<T>::addBondedTriMeshIG()\n";
951 console.
XDebug()<<
"BTriMeshIG name: " << name.c_str() <<
"\n";
953 console.
XDebug()<<
"BTriMeshIG mesh name: " << meshname.c_str() <<
"\n";
955 console.
XDebug()<<
"BTriMeshIG k : " << param.k <<
"\n";
957 console.
XDebug()<<
"BTriMeshIG r_break: " << param.brk <<
"\n";
959 console.
XDebug()<<
"BTriMeshIG build type: " << buildtype.c_str() <<
"\n";
963 if (m_mesh.find(meshname) != m_mesh.end())
965 tmp = m_mesh[meshname];
968 throw runtime_error(
"unknown mesh name in TSubLattice<T>::addTriMeshIG:" + meshname);
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"){
981 new_pis->buildFromPPAByGap(max_gap);
982 m_bpis.insert(make_pair(name,new_pis));
984 throw runtime_error(
"unknown build type in TSubLattice<T>::addBondedTriMeshIG:" + buildtype);
987 console.
XDebug() <<
"end TSubLattice<T>::addBondedTriMeshIG()\n";
996 console.
XDebug() <<
"TSubLattice<T>::addMesh2D()\n";
1001 vector<MeshNodeData2D> node_recv_buffer;
1002 vector<MeshEdgeData2D> edge_recv_buffer;
1009 console.
XDebug()<<
"Mesh2D name: " << name.c_str() <<
"\n";
1012 m_tml_comm.recv_broadcast_cont_packed(node_recv_buffer,0);
1013 console.
XDebug() <<
"recvd " << node_recv_buffer.size() <<
" nodes \n";
1016 m_tml_comm.recv_broadcast_cont_packed(edge_recv_buffer,0);
1017 console.
XDebug() <<
"recvd " << edge_recv_buffer.size() <<
" edges \n";
1021 new_tm->
LoadMesh(node_recv_buffer,edge_recv_buffer);
1023 m_mesh2d.insert(make_pair(name,new_tm));
1025 console.
XDebug() <<
"end TSubLattice<T>::addMesh2D()\n";
1035 console.
XDebug() <<
"TSubLattice<T>::addMesh2DIG()\n";
1045 console.
XDebug()<<
"Mesh2DIG type: " << type.c_str() <<
"\n";
1047 console.
XDebug()<<
"Mesh2DIG name: " << name.c_str() <<
"\n";
1049 console.
XDebug()<<
"Mesh2DIG mesh name: " << meshname.c_str() <<
"\n";
1053 if (m_mesh2d.find(meshname) != m_mesh2d.end())
1055 tmp = m_mesh2d[meshname];
1058 throw runtime_error(
"unknown mesh name in TSubLattice<T>::addMesh2DIG:" + meshname);
1067 m_dpis.insert(make_pair(name,new_pis));
1069 throw runtime_error(
"unknown type in TSubLattice<T>::addMesh2DIG:" + type);
1073 console.
XDebug() <<
"end TSubLattice<T>::addTriMeshIG()\n";
1083 console.
XDebug() <<
"TSubLattice<T>::addBondedMesh2DIG()\n";
1094 console.
XDebug() <<
"BMesh2DIG name: " << name.c_str() <<
"\n";
1096 console.
XDebug() <<
"BMesh2DIG mesh name: " << meshname.c_str() <<
"\n";
1098 console.
XDebug() <<
"BMesh2DIG k : " << param.k <<
"\n";
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";
1106 if (m_mesh2d.find(meshname) != m_mesh2d.end())
1108 tmp = m_mesh2d[meshname];
1111 throw runtime_error(
"unknown mesh name in TSubLattice<T>::addBondedMesh2DIG:" + meshname);
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"){
1125 m_bpis.insert(make_pair(name,new_pis));
1127 throw runtime_error(
"unknown build type in TSubLattice<T>::addBondedMesh2DIG:" + buildtype);
1130 console.
XDebug() <<
"end TSubLattice<T>::addBonded2DMeshIG()\n";
1142 console.
XDebug() <<
"TSubLattice<T>::addSingleIG()\n";
1149 console.
XDebug()<<
"SIG type: " << type.c_str() <<
"\n";
1152 if(type==
"Gravity"){
1156 m_singleParticleInteractions.insert(
1164 throw std::runtime_error(
1165 std::string(
"Trying to setup SIG of unknown type: ")
1171 console.
XDebug() <<
"end TSubLattice<T>::addSingleIG()\n";
1181 console.
XDebug() <<
"TSubLattice<T>::addDamping()\n";
1187 console.
XDebug()<<
"Damping type: " << type.c_str() <<
"\n";
1190 doAddDamping(type,param_buffer);
1192 console.
XDebug() <<
"end TSubLattice<T>::addDamping()\n";
1205 string damping_name;
1210 CDampingIGP *params=extractDampingIGP(¶m_buffer);
1212 damping_name=
"Damping";
1214 }
else if (type==
"ABCDamping"){
1217 damping_name=params->getName();
1219 }
else if (type==
"LocalDamping"){
1222 damping_name=params->getName();
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());
1234 m_damping.insert(make_pair(damping_name,DG));
1235 m_damping[damping_name]->update();
1247 console.
XDebug() <<
"TSubLattice<T>::addBondedIG()\n";
1253 param.tag=param_buffer.
pop_int();
1257 param.m_scaling=
static_cast<bool>(param_buffer.
pop_int());
1260 <<
"Got BondedIG parameters: " << param.tag
1261 <<
" " << name.c_str() <<
" "
1262 << param.
k <<
" " << param.
rbreak <<
"\n";
1270 console.
XDebug() <<
"set bpis unbreakable\"n";
1273 vector<int> vi(2,-1);
1274 for(
size_t i=0;i<m_temp_conn[param.tag].size();i+=2)
1276 vi[0] = (m_temp_conn[param.tag][i]);
1277 vi[1] = (m_temp_conn[param.tag][i+1]);
1278 B_PIS->tryInsert(vi);
1282 m_bpis.insert(make_pair(name,B_PIS));
1284 console.
XDebug() <<
"end TSubLattice<T>::addBondedIG()\n";
1294 console.
XDebug() <<
"TSubLattice<T>::addCappedBondedIG()\n";
1299 int tag=param_buffer.
pop_int();
1306 <<
"Got CappedBondedIG parameters: " << tag
1307 <<
" " << name.c_str() <<
" "
1308 << k <<
" " << rbreak <<
" " << maxforce <<
"\n";
1314 param.m_force_limit=maxforce;
1321 console.
XDebug() <<
"set bpis unbreakable\"n";
1335 vector<int> vi(2,-1);
1336 for(
size_t i=0;i<m_temp_conn[tag].size();i+=2)
1338 vi[0] = (m_temp_conn[tag][i]);
1339 vi[1] = (m_temp_conn[tag][i+1]);
1340 B_PIS->tryInsert(vi);
1344 m_bpis.insert(make_pair(name,B_PIS));
1346 console.
XDebug() <<
"end TSubLattice<T>::addCappedBondedIG()\n";
1352 console.
Error() <<
"TSubLattice<T>::addRotBondedIG() => trying to add rotational bonded IG to nonrotational model\n";
1358 console.
Error() <<
"TSubLattice<T>::addRotThermBondedIG() => trying to add rotational thermal bonded IG to nonrotational model\n";
1368 console.
XDebug() <<
"TSubLattice<T>::addShortBondedIG()\n";
1373 int tag=param_buffer.
pop_int();
1379 <<
"Got ShortBondedIG parameters: " << tag
1380 <<
" " << name.c_str() <<
" "
1381 << k <<
" " << rbreak <<
"\n";
1402 vector<int> vi(2,-1);
1403 for(
size_t i=0;i<m_temp_conn[param.tag].size();i+=2)
1405 vi[0] = (m_temp_conn[param.tag][i]);
1406 vi[1] = (m_temp_conn[param.tag][i+1]);
1407 B_PIS->tryInsert(vi);
1411 m_bpis.insert(make_pair(name,B_PIS));
1413 console.
XDebug() <<
"end TSubLattice<T>::addShortBondedIG()\n";
1424 console.
XDebug() <<
"TSubLattice<T>::addExIG()\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()))
1439 dynamic_ig->second->addExIG(bonded_ig->second);
1443 console.
Error() <<
"TSubLattice<T>::setExIG() - nonexisting interaction group \n";
1446 console.
XDebug() <<
"end TSubLattice<T>::addExIG()\n";
1457 console.
XDebug() <<
"TSubLattice<T>::removeIG()\n";
1463 map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.find(igname);
1464 if(iter!=m_dpis.end()){
1465 delete m_dpis[igname];
1468 console.
Error() <<
"TSubLattice<T>::removeIG() - nonexisting interaction group - ignore removal\n";
1479 console.
XDebug() <<
"TSubLattice<T>::exchangePos() \n" ;
1481 m_ppa->exchange(&T::getExchangeValues,&T::setExchangeValues);
1483 console.
XDebug() <<
"end TSubLattice<T>::exchangePos()\n" ;
1492 console.
XDebug() <<
"TSubLattice<T>::zeroForces()\n";
1495 m_ppa->forAllParticles(&T::zeroForce);
1497 for(map<string,TriMesh*>::iterator iter=m_mesh.begin();
1500 (iter->second)->zeroForces();
1504 for(map<string,Mesh2D*>::iterator iter=m_mesh2d.begin();
1505 iter!=m_mesh2d.end();
1507 (iter->second)->zeroForces();
1512 (iter->second)->zeroForce();
1514 console.
XDebug() <<
"end TSubLattice<T>::zeroForces() \n";
1526 console.
XDebug() <<
"TSubLattice<T>::calcForces() \n";
1531 (it->second)->calcForces();
1535 typename NameIGroupMap::iterator siter=this->m_singleParticleInteractions.begin();
1536 siter != m_singleParticleInteractions.end();
1540 (siter->second)->calcForces();
1543 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();iter!=m_dpis.end();iter++)
1545 (iter->second)->calcForces();
1548 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();iter!=m_bpis.end();iter++)
1550 (iter->second)->calcForces();
1553 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();iter!=m_damping.end();iter++)
1555 (iter->second)->calcForces();
1558 console.
XDebug() <<
"end TSubLattice<T>::calcForces() \n";
1571 console.
XDebug() <<
"TSubLattice<T>::setTimeStepSize() \n";
1576 (it->second)->setTimeStepSize(dt);
1580 typename NameIGroupMap::iterator siter=this->m_singleParticleInteractions.begin();
1581 siter != m_singleParticleInteractions.end();
1585 (siter->second)->setTimeStepSize(dt);
1588 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();iter!=m_dpis.end();iter++)
1590 (iter->second)->setTimeStepSize(dt);
1593 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();iter!=m_bpis.end();iter++)
1595 (iter->second)->setTimeStepSize(dt);
1598 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();iter!=m_damping.end();iter++)
1600 (iter->second)->setTimeStepSize(dt);
1603 console.
XDebug() <<
"end TSubLattice<T>::setTimeStepSize() \n";
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";
1630 if (this->getParticleType() ==
"RotTherm")
1632 this->oneStepTherm();
1645 integrateTherm(m_dt);
1657 console.
XDebug() <<
"TSubLattice<T>::integrateTherm \n";
1658 m_ppa->forAllParticles(&T::integrateTherm,dt);
1660 console.
XDebug() <<
"end TSubLattice<T>::integrateTherm \n";
1666 console.
XDebug() <<
"TSubLattice<T>::thermExpansion() \n";
1667 m_ppa->forAllParticles(&T::thermExpansion);
1669 console.
XDebug() <<
"end TSubLattice<T>::thermExpansion() \n";
1678 console.
XDebug() <<
"TSubLattice<T>::zeroHeat()\n";
1681 m_ppa->forAllParticles(&T::zeroHeat);
1690 console.
XDebug() <<
"end TSubLattice<T>::zeroHeat() \n";
1702 console.
XDebug() <<
"TSubLattice<T>::calcHeatFrict() \n";
1706 typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1711 (iter->second)->calcHeatFrict();
1714 console.
XDebug() <<
"end TSubLattice<T>::calcHeatFrict() \n";
1720 console.
XDebug() <<
"TSubLattice<T>::calcHeatTrans() \n";
1725 typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1730 (iter->second)->calcHeatTrans();
1734 typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();
1739 (iter->second)->calcHeatTrans();
1742 console.
XDebug() <<
"end TSubLattice<T>::calcHeatTrans() \n";
1761 m_pTimers->start(
"RebuildInteractions");
1762 m_pTimers->resume(
"NeighbourSearch");
1763 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();
1767 console.
Debug() <<
"exchg & rebuild BPIS " << iter->first <<
" at node " << m_rank <<
"\n";
1768 (iter->second)->exchange();
1769 (iter->second)->rebuild();
1772 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
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();
1785 resetDisplacements();
1786 m_pTimers->stop(
"RebuildInteractions");
1795 console.
Debug() <<
"CSubLattice<T>::searchNeighbors()\n";
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";
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;
1820 m_pTimers->start(
"UpdateBondedInteractions");
1821 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_bpis.begin();
1825 bool n=(iter->second)->update();
1826 need_update=need_update || n;
1828 m_pTimers->stop(
"UpdateBondedInteractions");
1829 if((m_ppa->getTimeStamp() > m_last_ns) || need_update)
1831 m_pTimers->start(
"UpdateDynamicInteractions");
1832 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
1836 bool n=(iter->second)->update();
1837 need_update=need_update || n;
1839 m_pTimers->stop(
"UpdateDynamicInteractions");
1844 (it->second)->Update(m_ppa);
1846 for(
typename map<string,AParallelInteractionStorage*>::iterator iter=m_damping.begin();
1847 iter!=m_damping.end();
1849 (iter->second)->update();
1851 m_last_ns=m_ppa->getTimeStamp();
1854 console.
Debug() <<
"end TSubLattice<T>::updateInteractions()\n";
1864 console.
Debug() <<
"TSubLattice<T>::checkNeighbors()\n";
1867 double alpha=0.5*m_alpha;
1868 double srsqr=alpha*alpha;
1874 m_ppa->forAllParticlesGet(displ,&T::getDisplacement);
1877 vector<Vec3>::iterator it=displ.begin();
1878 while((it!=displ.end())&&(mdsqr<srsqr))
1880 double sqdisp=(*it)*(*it);
1881 mdsqr = ((mdsqr < sqdisp) ? sqdisp : mdsqr);
1885 console.
XDebug() <<
"max squared displacement " << mdsqr <<
"alpha^2 = " << srsqr <<
"\n";
1895 for(map<string,TriMesh*>::iterator iter=m_mesh.begin();
1899 if(iter->second->hasMovedBy(alpha)){
1906 console.
Debug() <<
"end TSubLattice<T>::checkNeighbors()\n";
1916 console.
Debug() <<
"slave " << m_rank <<
" resetDisplacements()\n";
1917 m_ppa->forAllParticles(&T::resetDisplacement);
1918 for(map<string,TriMesh*>::iterator iter=m_mesh.begin();
1921 iter->second->resetCurrentDisplacement();
1923 console.
Debug() <<
"slave " << m_rank <<
" end resetDisplacements()\n";
1933 console.
Debug() <<
"TSubLattice<T>::moveParticleTo()\n";
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";
1950 console.
Debug() <<
"TSubLattice<T>::moveTaggedParticlesBy()\n";
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";
1964 m_ppa->forParticle(particleId, (
void (T::*)(
Vec3))(&T::moveTo), posn);
1974 console.
Debug() <<
"TSubLattice<T>::moveSingleNode()\n";
1980 Vec3 disp=buffer.pop_vector();
1982 console.
XDebug() <<
"name :" << name <<
" id : " <<
id <<
" disp " << disp <<
"\n";
1984 map<string,TriMesh*>::iterator tm=m_mesh.find(name);
1985 if (tm!=m_mesh.end()){
1986 (tm->second)->moveNode(
id,disp);
1988 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(name);
1989 if(m2d!=m_mesh2d.end()){
1990 (m2d->second)->moveNode(
id,disp);
1993 console.
Debug() <<
"end TSubLattice<T>::moveSingleNode()\n";
2003 console.
Error() <<
"TSubLattice<T>::moveTaggedNodes() NOT IMPLEMENTED\n";
2006 "TSubLattice<T>::moveTaggedNodes() NOT IMPLEMENTED\n"
2014 Vec3 disp=buffer.pop_vector();
2027 map<string,TriMesh*>::iterator tm=m_mesh.find(meshName);
2028 if (tm != m_mesh.end()){
2029 (tm->second)->translateBy(translation);
2031 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshName);
2032 if(m2d!=m_mesh2d.end()){
2033 (m2d->second)->translateBy(translation);
2040 console.
Debug() <<
"TSubLattice<T>::findParticleNearestTo: enter\n";
2041 const T *pClosest = NULL;
2042 double minDistSqrd = std::numeric_limits<double>::max();
2044 typename ParticleArray::ParticleIterator it =
2045 m_ppa->getInnerParticleIterator();
2046 while (it.hasNext())
2048 const T &p = it.next();
2049 const double distSqrd = (pt - p.getPos()).norm2();
2050 if (distSqrd < minDistSqrd)
2052 minDistSqrd = distSqrd;
2056 console.
Debug() <<
"TSubLattice<T>::findParticleNearestTo: exit\n";
2061 std::make_pair(sqrt(minDistSqrd), pClosest->getID())
2063 std::make_pair(std::numeric_limits<double>::max(), -1)
2073 const T *particle = NULL;
2075 m_ppa->getInnerParticleIterator();
2076 while (it.hasNext())
2078 const T &p = it.next();
2079 if (p.getID() == particleId)
2084 if (particle != NULL)
2086 return std::make_pair(particleId, particle->getPos());
2088 return std::make_pair(-1,Vec3::ZERO);
2098 <<
"TSubLattice<T>::getParticleData: enter\n";
2099 typedef std::set<int> IdSet;
2100 typedef std::vector<T> ParticleVector;
2102 ParticleVector particleVector;
2104 m_ppa->getInnerParticleIterator();
2105 if (particleIdVector.size() > 0)
2107 IdSet idSet(particleIdVector.begin(), particleIdVector.end());
2109 <<
"TSubLattice<T>::getParticleData: iterating over particles\n";
2110 while (it.hasNext())
2112 const T &p = it.next();
2113 if (idSet.find(p.getID()) != idSet.end())
2115 particleVector.push_back(p);
2121 m_ppa->getAllInnerParticles(particleVector);
2124 <<
"TSubLattice<T>::getParticleData:"
2125 <<
" sending particle data of size " << particleVector.size() <<
"\n";
2126 m_tml_comm.send_gather_packed(particleVector, 0);
2128 <<
"TSubLattice<T>::getParticleData: exit\n";
2137 console.
Debug() <<
"TSubLattice<T>::tagParticleNearestTo()\n";
2143 Vec3 pos=buffer.pop_vector();
2146 T* part_ptr=m_ppa->getParticlePtrByPosition(pos);
2148 int old_tag=part_ptr->getTag();
2149 int new_tag=(old_tag & (~mask)) | (tag & mask);
2150 part_ptr->setTag(new_tag);
2152 cout <<
"pos, realpos: " << pos <<
" " << part_ptr->getPos() <<
" old tag, new tag " << old_tag <<
" " << part_ptr->getTag() << endl;
2154 console.
Debug() <<
"end TSubLattice<T>::tagParticleNearestTo()\n";
2164 console.
Debug() <<
"TSubLattice<T>::setParticleNonDynamic()\n";
2169 m_ppa->forParticleTag(tag,(
void (T::*)())(&T::setNonDynamic));
2170 console.
Debug() <<
"end TSubLattice<T>::setParticleNonDynamic()\n";
2180 console.
Debug() <<
"TSubLattice<T>::setParticleNonRot()\n";
2185 m_ppa->forParticleTag(tag,(
void (T::*)())(&T::setNonDynamicRot));
2186 console.
Debug() <<
"end TSubLattice<T>::setParticleNonRot()\n";
2196 console.
Debug() <<
"TSubLattice<T>::setParticleNonTrans()\n";
2201 m_ppa->forParticleTag(tag,(
void (T::*)())(&T::setNonDynamicLinear));
2202 console.
Debug() <<
"end TSubLattice<T>::setParticleNonTrans()\n";
2211 console.
Debug() <<
"TSubLattice<T>::setTaggedParticleVel()\n";
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";
2227 console.
XDebug() <<
"TSubLattice<T>::moveWallBy()\n";
2232 Vec3 mv=buffer.pop_vector();
2233 typename map<string,CWall*>::iterator iter=m_walls.find(wname);
2234 if(iter!=m_walls.end())
2236 (iter->second)->moveBy(mv);
2246 console.
XDebug() <<
"TSubLattice<T>::setWallNormal()\n";
2251 Vec3 wn=buffer.pop_vector();
2252 typename map<string,CWall*>::iterator iter=m_walls.find(wname);
2253 if(iter!=m_walls.end())
2255 (iter->second)->setNormal(wn);
2265 console.
XDebug() <<
"TSubLattice<T>::applyForceToWall()\n";
2270 Vec3 f=buffer.pop_vector();
2271 typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.find(wname);
2272 if(iter!=m_WIG.end())
2274 (iter->second)->applyForce(f);
2285 console.
XDebug() <<
"TSubLattice<T>::setVelocityOfWall()\n";
2290 Vec3 v=buffer.pop_vector();
2291 typename map<string,AWallInteractionGroup<T>*>::iterator iter=m_WIG.find(wname);
2292 if(iter!=m_WIG.end())
2294 (iter->second)->setVelocity(v);
2304 console.
Debug() <<
"TSubLattice<T>::setParticleVelocity()\n";
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";
2320 console.
Debug() <<
"TSubLattice<T>::setParticleDensity()\n";
2327 m_ppa->forParticleTagMask(tag,mask,(
void (T::*)(
double))(&T::setDensity),rho);
2328 console.
XDebug() <<
"end TSubLattice<T>::setParticleVelocity()\n";
2338 console.
Debug() <<
"TSubLattice<T>::sendDataToMaster()\n";
2339 vector<Vec3> positions;
2340 vector<double> radii;
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));
2347 m_tml_comm.send_gather(positions,0);
2348 m_tml_comm.send_gather(radii,0);
2349 m_tml_comm.send_gather(ids,0);
2351 console.
Debug() <<
"end TSubLattice<T>::sendDataToMaster()\n";
2360 console.
Debug()<<
"TSubLattice<T>::countParticles()\n";
2363 buffer.
append(m_ppa->size());
2374 cout<<
"My Rank : " << m_rank <<
"\n" ;
2377 cout << *m_ppa << endl;
2384 cout <<
"Data: my rank : " << m_rank <<
"particles : \n" ;
2385 m_ppa->forAllParticles((
void (T::*)())(&T::print));
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";
2411 m_tml_comm.recv_broadcast_cont(fieldname,0);
2413 m_tml_comm.recv_broadcast(
id,0);
2415 m_tml_comm.recv_broadcast(is_tagged,0);
2418 typename T::ScalarFieldFunction rdf=T::getScalarFieldFunction(fieldname);
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";
2433 m_field_slaves.insert(make_pair(
id,new_spfs));
2442 console.
XDebug() <<
"TSubLattice<T>::addVectorParticleField\n";
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";
2453 typename T::VectorFieldFunction rdf=T::getVectorFieldFunction(fieldname);
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";
2468 m_field_slaves.insert(make_pair(
id,new_vpfs));
2470 console.
Debug() <<
"end TSubLattice<T>::addVectorParticleField\n";
2480 console.
XDebug() <<
"TSubLattice<T>::addScalarInteractionField\n";
2484 int id,is_tagged,tag,mask,is_checked;
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";
2498 map<string,AParallelInteractionStorage*>::iterator it=m_dpis.find(igname);
2501 m_tml_comm.recv_broadcast(tag,0);
2502 m_tml_comm.recv_broadcast(mask,0);
2504 m_tml_comm.recv_broadcast(is_checked,0);
2505 console.
XDebug() <<
"recvd. is_checked: " << is_checked <<
"\n";
2507 if(it!=m_dpis.end())
2509 console.
XDebug() <<
"found interaction group in dynamic\n";
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));
2516 it=m_bpis.find(igname);
2517 if(it!=m_bpis.end()){
2518 console.
XDebug() <<
"found interaction group in bonded\n";
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));
2526 it=m_damping.find(igname);
2527 if(it!=m_damping.end())
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));
2535 cerr <<
"ERROR : unknown interaction group name in addScalarInteractionField " << endl;
2540 console.
XDebug() <<
"end TSubLattice<T>::addScalarInteractionField\n";
2549 console.
Debug() <<
"TSubLattice<T>::addVectorInteractionField\n";
2553 int id,is_tagged,tag,mask,is_checked;
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";
2567 map<string,AParallelInteractionStorage*>::iterator it=m_dpis.find(igname);
2570 m_tml_comm.recv_broadcast(tag,0);
2571 m_tml_comm.recv_broadcast(mask,0);
2573 m_tml_comm.recv_broadcast(is_checked,0);
2574 console.
XDebug() <<
"recvd. is_checked: " << is_checked <<
"\n";
2576 if(it!=m_dpis.end())
2578 console.
XDebug() <<
"found interaction group in dynamic\n";
2580 new_sifs=(it->second)->generateNewVectorFieldSlave(&m_tml_comm,fieldname,is_checked,is_tagged,tag,mask);
2582 m_field_slaves.insert(make_pair(
id,new_sifs));
2584 console.
Error()<<
"ERROR: could not generate Field Slave for field " << fieldname <<
"\n";
2589 it=m_bpis.find(igname);
2590 if(it!=m_bpis.end()){
2591 console.
XDebug() <<
"found interaction group in bonded\n";
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));
2599 it=m_damping.find(igname);
2600 if(it!=m_damping.end())
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));
2608 cerr <<
"ERROR : unknown interaction group name in addScalarInteractionField " << endl;
2613 console.
Debug() <<
"end TSubLattice<T>::addVectorInteractionField\n";
2623 console.
Debug() <<
"TSubLattice<T>::addVectorTriangleField()\n";
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";
2636 map<string,TriMesh*>::iterator tm=m_mesh.find(meshname);
2638 if (tm!=m_mesh.end()){
2643 console.
Critical() <<
"NULL rdf !!!\n";
2646 m_field_slaves.insert(make_pair(
id,new_vfs));
2648 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname);
2649 if(m2d!=m_mesh2d.end()){
2653 console.
Critical() <<
"NULL rdf !!!\n";
2656 m_field_slaves.insert(make_pair(
id,new_efs));
2659 console.
Debug() <<
"end TSubLattice<T>::addVectorTriangleField()\n";
2668 console.
Debug() <<
"TSubLattice<T>::addScalarTriangleField()\n";
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";
2685 console.
Critical() <<
"NULL rdf !!!\n";
2688 m_field_slaves.insert(make_pair(
id,new_vtfs));
2689 console.
Debug() <<
"end TSubLattice<T>::addScalarTriangleField()\n";
2698 console.
XDebug() <<
"begin TSubLattice<T>::addVectorWallField()\n";
2700 string tmp_wallname;
2701 vector<string> wallnames;
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();
2716 m_tml_comm.recv_broadcast(
id,0);
2717 console.
XDebug() <<
"recvd. id: " <<
id <<
"\n";
2720 map<string,CWall*>::iterator cwalliter=m_walls.find(*(wallnames.begin()));
2721 if(cwalliter==m_walls.end()){
2722 std::stringstream 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++)
2728 msg <<
"'" << it->first <<
"' ";
2730 console.
Error() << msg.str() <<
"\n";
2731 throw std::runtime_error(msg.str());
2734 int sumflag=(cwalliter->second)->getFieldSummationFlag(fieldname);
2736 if(m_tml_comm.rank()==1){
2737 m_tml_comm.send(sumflag,0);
2739 m_tml_comm.barrier();
2742 AWallFieldSlave* new_fs=(cwalliter->second)->generateVectorFieldSlave(&m_tml_comm,fieldname);
2745 vector<string>::iterator niter=wallnames.begin();
2746 if(niter!=wallnames.end()) niter++ ;
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()){
2751 std::stringstream 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++)
2757 msg <<
"'" << it->first <<
"' ";
2760 console.
Error() << msg.str() <<
"\n";
2761 throw std::runtime_error(msg.str());
2763 new_fs->
addWall(iter->second);
2768 m_field_slaves.insert(make_pair(
id,new_fs));
2770 console.
Error() <<
"ERROR in addVectorWallField: got NULL Slave\n";
2774 console.
XDebug() <<
"end TSubLattice<T>::addVectorWallField()\n";
2783 console.
Debug() <<
"TSubLattice<T>::sendFieldData()\n";
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)
2790 m_field_slaves[id]->sendData();
2794 cerr <<
"NULL pointer in m_field_slaves!" << endl;
2797 console.
Debug() <<
"end TSubLattice<T>::sendFieldData()\n";
2809 std::streamsize oldprec=oStream.precision(9);
2816 particleIt(particleArray.getInnerParticleIterator());
2818 const std::string delim =
"\n";
2820 oStream << particleIt.getNumRemaining() << delim;
2821 while (particleIt.hasNext()) {
2822 particleIt.next().saveSnapShotData(oStream);
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);
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();
2842 oStream << tm_iter->first << delim;
2843 tm_iter->second->writeCheckPoint(oStream,delim);
2847 oStream.precision(oldprec);
2857 std::streamsize oldprec=oStream.precision(16);
2859 const std::string delim =
"\n";
2864 m_ppa->saveCheckPointData(oStream);
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);
2881 for(std::map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
2884 if(iter->second->willSave()) count_save++;
2886 oStream << count_save << delim;
2887 for(std::map<string,AParallelInteractionStorage*>::iterator iter=m_dpis.begin();
2890 if(iter->second->willSave()) {
2891 iter->second->saveCheckPointData(oStream);
2897 oStream <<
"Walls " << m_walls.size() << delim;
2898 for(map<string,CWall*>::iterator w_iter=m_walls.begin();
2899 w_iter!=m_walls.end();
2901 oStream << w_iter->first << delim;
2902 w_iter->second->writeCheckPoint(oStream,delim);
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();
2910 oStream << tm_iter->first << delim;
2911 tm_iter->second->writeCheckPoint(oStream,delim);
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();
2918 oStream << tm_iter->first << delim;
2919 tm_iter->second->writeCheckPoint(oStream,delim);
2922 oStream.precision(oldprec);
2929 m_ppa->loadCheckPointData(iStream);
2934 barrier.wait(
"PPA rebuild");
2939 iStream >> nr_bonded_ig;
2943 if(nr_bonded_ig!=m_bpis.size()){
2944 std::cerr <<
"number of bonded interaction groups differ between snapshot and script!" << std::endl;
2946 for (map<string,AParallelInteractionStorage*>::iterator it = m_bpis.begin();
2949 it->second->loadCheckPointData(iStream);
2954 int nr_nonbonded_ig;
2955 iStream >> nr_nonbonded_ig;
2959 if(nr_nonbonded_ig!=m_dpis.size()){
2960 std::cerr <<
"number of dynamic interaction groups differ between snapshot and script!" << std::endl;
2962 for (map<string,AParallelInteractionStorage*>::iterator it = m_dpis.begin();
2965 it->second->loadCheckPointData(iStream);
2971 if(token!=
"Walls") {
2972 std::cerr <<
"expected Walls , got " << token << std::endl;
2979 for(
int i=0;i<nwalls;i++){
2983 m_walls[wname]=new_wall;
2990 if(token!=
"TriMesh") {
2991 std::cerr <<
"expected TriMesh , got " << token << std::endl;
2996 for(
int i=0;i<nmesh;i++){
2999 new_tm->loadCheckPoint(iStream);
3000 m_mesh.insert(make_pair(mname,new_tm));
3004 if(token!=
"Mesh2D") {
3005 std::cerr <<
"expected Mesh2D , got " << token << std::endl;
3010 for(
int i=0;i<nmesh;i++){
3014 m_mesh2d.insert(make_pair(mname,new_m2d));
3026 console.
XDebug() <<
"TSubLattice<T>::getMeshNodeRef()\n";
3027 vector<int> ref_vec;
3034 console.
XDebug() <<
"Mesh name: " << meshname <<
"\n";
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();
3042 ref_vec.push_back(iter->getID());
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();
3050 ref_vec.push_back(iter->getID());
3053 console.
Critical() <<
"ERROR - WRONG MESH NAME IN getMeshNodeRef() !! \n";
3057 m_tml_comm.send_gather(ref_vec,0);
3059 console.
XDebug() <<
"end TSubLattice<T>::getMeshNodeRef()\n";
3068 console.
XDebug() <<
"TSubLattice<T>::getMeshFaceRef()\n";
3069 vector<int> ref_vec;
3076 console.
XDebug() <<
"Mesh name: " << meshname <<
"\n";
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();
3084 ref_vec.push_back(iter->getID());
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();
3092 ref_vec.push_back(iter->getID());
3095 console.
Critical() <<
"ERROR - WRONG MESH NAME IN getMeshNodeRef() !! \n";
3099 m_tml_comm.send_gather(ref_vec,0);
3101 console.
XDebug() <<
"end TSubLattice<T>::getMeshNodeRef()\n";
3110 console.
XDebug() <<
"TSubLattice<T>::getMesh2DStress()\n";
3116 console.
XDebug() <<
"Mesh name: " << meshname <<
"\n";
3119 map<string,Mesh2D*>::iterator m2d=m_mesh2d.find(meshname);
3120 if(m2d!=m_mesh2d.end()){
3121 vector<pair<int,Vec3> > data_vec;
3123 data_vec=(m2d->second)->forAllEdgesGetIndexed(&Edge2D::getForceDensity);
3126 m_tml_comm.send_gather(data_vec,0);
3128 console.
Critical() <<
"ERROR - WRONG MESH NAME IN getMesh2DStress() !! \n";
3131 console.
XDebug() <<
"end TSubLattice<T>::getMesh2DStress()\n";
3140 console.
XDebug() <<
"TSubLattice<T>::getTriMeshStress(): enter\n";
3145 const std::string meshName = param_buffer.
pop_string();
3146 console.
XDebug() <<
"Mesh name: " << meshName <<
"\n";
3149 map<string,TriMesh*>::iterator m=m_mesh.find(meshName);
3150 if(m != m_mesh.end()){
3151 vector<pair<int,Vec3> > data_vec;
3153 data_vec=(m->second)->forAllTrianglesGetIndexed(&Triangle::getForce);
3156 m_tml_comm.send_gather(data_vec,0);
3158 std::stringstream msg;
3159 msg <<
"Invalid mesh name: " << meshName <<
". No such triangular mesh.";
3160 throw std::runtime_error(msg.str().c_str());
3163 console.
XDebug() <<
"TSubLattice<T>::getTriMeshStress(): exit\n";