Allink  v0.1
ForcesBoundary.cpp
1 #include "Forces.h"
2 
4  Pm[Bead2Move].Pos[2] += IncrDist;
5  if(VAR_IF_TYPE(SysShape,SYS_ROD))
6  Pm[Bead2Move-1].Pos[2] += IncrDist;
7 }
9  Pm[Bead2Move].Pos[2] -= IncrDist;
10  if(VAR_IF_TYPE(SysShape,SYS_ROD))
11  Pm[Bead2Move-1].Pos[2] -= IncrDist;
12 }
13 void Forces::SelectBead(int p){
14  if(p<0 && p > pNPart()) return;
15  Bead2Move = p;
16  if(pType(p) != 2)
17  Pm[p].Typ = 1;
18  if(pType(Old2Move) == 1)
19  Pm[Old2Move].Typ = 0;
20  Old2Move = p;
21 }
22 void Forces::Wave(){
23  for(int s=0;s<NEdge;s++){
24  for(int ss=0;ss<NEdge;ss++){
25  Pm[s*NEdge+ss].Pos[2] = .5+.1*sin(s*Dx*60.*INV_DUE_PI+Time)*sin(ss*Dx*30.*INV_DUE_PI+Time);
26  }
27  }
28 }
29 void Forces::AddCircle(int n){
30  for(int p=0;p<pNPart();p++){
31  double CosTh = (Nano[n].Pos[0] - pPos(p,0))/Nano[n].Rad;
32  double SenTh = sqrt(1 - QUAD(CosTh));
33  double Sign = 1.;
34  for(int d=0;d<3;d++){
35  Pm[p].Vel[d] += (Fm[p].Dir[d] + Fm[p].Ext[d])*pDeltat();
36  }
37  if(pType(p) != 2) Pm[p].Typ = 0;
38  if(p <= NEdge){
39  Sign *= -1.;
40  if(Pm[p].Pos[2] + .5*Pm[p].Vel[2]*pDeltat() >= Sign*(Nano[n].Rad)*SenTh + Nano[n].Pos[2]){
41  Pm[p].Pos[2] = Sign*(Nano[n].Rad)*SenTh + Nano[n].Pos[2];
42  Pm[p].Vel[2] = 0.;
43  Sign = 1.;
44  if(Pm[p].Pos[0] > Nano[n].Pos[0])
45  Sign = -1.;
46  if(p > 1)
47  Fm[p-1].Dir[0] -= Kf.Cont*SenTh*Sign;
48  if(p < NEdge - 1)
49  Fm[p+1].Dir[0] += Kf.Cont*SenTh*Sign;
50  if(Pm[p].Typ != 2)
51  Pm[p].Typ = 1;
52  }
53  else
54  Fm[p].Dir[0] = 0.;
55  }
56  else if( p> NEdge){
57  if(Pm[p].Pos[2] + .5*Pm[p].Vel[2]*pDeltat() <= Sign*(Nano[n].Rad)*SenTh + Nano[n].Pos[2]){
58  Pm[p].Pos[2] = Sign*(Nano[n].Rad)*SenTh + Nano[n].Pos[2];
59  Pm[p].Vel[2] = 0.;
60  Sign = 1.;
61  if(Pm[p].Pos[0] > Nano[n].Pos[0])
62  Sign = -1.;
63  if(p > NEdge)
64  Fm[p-1].Dir[0] -= Kf.Cont*SenTh*Sign;
65  if(p < pNPart() - 1)
66  Fm[p+1].Dir[0] += Kf.Cont*SenTh*Sign;
67  if(Pm[p].Typ != 2)
68  Pm[p].Typ = 1;
69  }
70  else
71  Fm[p].Dir[0] = 0.;
72  }
73  }
74 }
75 void Forces::AddCylinder(int n){
76  for(int p=0;p<pNPart();p++){
77  if(Pm[p].Typ == 2) continue;
78  double Dist = sqrt(SQR(Pm[p].Pos[0]-Nano[n].Pos[0])+SQR(Pm[p].Pos[1]-Nano[n].Pos[1]));
79  if(Dist > Nano[n].Rad) continue;
80  double Zed = (Nano[n].Rad-Dist)*tan(Nano[n].Hamaker*DUE_PI/360.);
81  Pm[p].Typ = 1;
82  Pm[p].Pos[2] = .5*Nano[n].Height + Nano[n].Pos[2] + Zed;
83  }
84 }
88 double Forces::HeightBoundary(double *Pos,int dir){
89  return Pos[2];
90  int nNano = 0;
91  int dir2 = (dir+1)%2;
92  for(int n=0;n<pNNano();n++){
93  if(fabs(Pos[dir] - Nano[n].Pos[dir]) < Nano[n].Rad){
94  nNano = n;
95  break;
96  }
97  }
98  double yNano = Pos[dir2] - Nano[nNano].Pos[dir2];
99  double Lat1 = Nano[nNano].Pos[dir] - sqrt(SQR(Nano[nNano].Rad)-SQR(yNano) ) - Pos[dir];
100  double Lat2 = Nano[nNano].Pos[dir] + sqrt(SQR(Nano[nNano].Rad)-SQR(yNano) ) - Pos[dir];
101  double Dist = MIN(fabs(Lat1),fabs(Lat2));
102  //double Dist = sqrt(SQR(Pos[0]-Nano[nNano].Pos[0])+SQR(Pos[1]-Nano[nNano].Pos[1]));
103  //double Zed = (Nano[nNano].Rad-Dist)*tan(Nano[nNano].Hamaker*DUE_PI/360.);
104  double Zed = Dist*tan(Nano[nNano].Hamaker*DUE_PI/360.);
105  Pos[2] = .5*Nano[nNano].Height + Nano[nNano].Pos[2] + Zed;
106  return .5*Nano[nNano].Height + Nano[nNano].Pos[2] + Zed;
107 }
108 void Forces::AddPore(int n){
109  for(int p=0;p<pNPart();p++){
110  if(Pm[p].Typ == 2) continue;
111  double Dist = sqrt(SQR(Pm[p].Pos[0]-Nano[n].Pos[0])+SQR(Pm[p].Pos[1]-Nano[n].Pos[1]));
112  if(Dist > Nano[n].Rad) continue;
113  Pm[p].Typ = 1;
114  Pm[p].Pos[2] = Nano->Pos[2];
115  }
116 }
118  for(int p=0;p<pNPart();p++){
119  if(Pm[p].Typ == 2) continue;
120  Pm[p].Typ = 0;
121  }
122  for(int n=0;n<pNNano();n++){
123  if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_SPH)){
124  AddCircle(n);
125  }
126  else if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_CYL)){
127  AddCylinder(n);
128  }
129  else if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_PORE)){
130  AddPore(n);
131  }
132  }
133 }
double Cont
Contact/Friction.
Definition: Forces.h:157
double Vel[4]
xyzr Velocity of the particle
Definition: VarData.h:220
void AddRigid()
Add all rigid bodies as a boundary condition.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
void PullBead()
Move a particle.
double Rad
Size.
Definition: VarData.h:445
double HeightBoundary(double *Pos, int dir)
Height of the boundary condition depending on the direction.
double Height
Height of the cylinder.
Definition: VarData.h:449
int pType(int p)
Return the type.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
double Dir[3]
Direction.
Definition: Forces.h:178
double Pos[3]
Position.
Definition: VarData.h:427
double Time
Total time.
Definition: Forces.h:733
double Dx
Spatial separation between particles.
Definition: Forces.h:729
void PushBead()
Move a particle.
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
double IncrDist
Step to move a particle.
Definition: Forces.h:715
int Old2Move
Old part to move.
Definition: Forces.h:739
int SysShape
Shape of system.
Definition: Forces.h:745
int pNNano()
Number of nanoparticles.
void AddPore(int nNano)
Add a pore as a boundary condition.
void Wave()
Sinusoidal surface wave.
double pPos(int p, int d)
Return back folded position.
void AddCylinder(int nNano)
Add a cylinder as a boundary condition.
double Ext[3]
External.
Definition: Forces.h:182
void SelectBead(int p)
Select a particle.
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int Typ
Type.
Definition: VarData.h:226
double pDeltat()
Delta t.
Definition: VarData.h:910
void AddCircle(int nNano)
Add a circle as a boundary condition.
FORCES * Fm
Array containing the forces for each particle.
Definition: Forces.h:781
int Bead2Move
Bead to move.
Definition: Forces.h:737
int pNPart()
Number of particle.