1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
/////////////////////////////////////////////////////////////
// //
// Copyright (c) 2003-2014 by The University of Queensland //
// Centre for Geoscience Computing //
// http://earth.uq.edu.au/centre-geoscience-computing //
// //
// Primary Business: Brisbane, Queensland, Australia //
// Licensed under the Open Software License version 3.0 //
// http://www.opensource.org/licenses/osl-3.0.php //
// //
/////////////////////////////////////////////////////////////
#include "LinearDashpotInteraction.h"
// --- system includes ---
#include <iostream>
// --- member functions for interaction group parameter class ---
//! default constructor
CLinearDashpotIGP::CLinearDashpotIGP(): AIGParam(), m_damp(0.0), m_cutoff(1.0)
{}
/*!
constructor with parameters
\param name the name of the interaction group
\param damp the damping coefficient
\param cutoff the interaction range, relative to the sum of the particle radii
*/
CLinearDashpotIGP::CLinearDashpotIGP(const std::string& name,double damp, double cutoff)
: AIGParam(name), m_damp(damp), m_cutoff(cutoff)
{}
// --- function of interaction class ---
CLinearDashpotInteraction::CLinearDashpotInteraction(CParticle* p1,CParticle* p2, const CLinearDashpotIGP& param) : APairInteraction(p1,p2)
{
// calc. cross section - 2D / 3D
double r_avg=0.5*(m_p1->getRad()+m_p2->getRad());
if(CParticle::getDo2dCalculations()){
m_cross_section=2.0*r_avg;
} else {
m_cross_section=M_PI*r_avg*r_avg;
}
// set parameters
m_cutoff=param.m_cutoff;
m_damp=param.m_damp;
m_force=Vec3(0.0,0.0,0.0);
}
/*!
calculate forces
*/
void CLinearDashpotInteraction::calcForces()
{
// calc distance -> needed to check if particles interact
Vec3 D=m_p1->getPos()-m_p2->getPos();
double dist_sq=D*D;
// calc cutoff distance -> could be cached
double cut_dist=m_cutoff*(m_p1->getRad()+m_p2->getRad());
if(dist_sq<(cut_dist*cut_dist)){
// velocity difference
Vec3 dvel=m_p1->getVel()-m_p2->getVel();
// strain rate
Vec3 eps_dot=dvel/sqrt(dist_sq);
m_force=eps_dot*m_damp*m_cross_section;
// apply force at particle centers
m_p2->applyForce(m_force,m_p2->getPos());
m_p1->applyForce(-1.0*m_force,m_p1->getPos());
}
m_cpos=(m_p1->getPos()+m_p2->getPos())*0.5;
}
/*!
"field function" returning force currently exerted by interaction
*/
Vec3 CLinearDashpotInteraction::getForce() const
{
return m_force;
}
/*!
Get the particle member function which returns a scalar field of a given name.
\param name the name of the field
*/
CLinearDashpotInteraction::ScalarFieldFunction CLinearDashpotInteraction::getScalarFieldFunction(const string& name)
{
CLinearDashpotInteraction::ScalarFieldFunction sf;
if (name=="count"){
sf=&CLinearDashpotInteraction::Count;
} else {
sf=NULL;
std::cerr << "ERROR - invalid name for interaction scalar access function " << name << " in LinearDashpotInteraction" << std::endl;
}
return sf;
}
/*!
Get the particle member function which returns a vector field of a given name.
\param name the name of the field
*/
CLinearDashpotInteraction::VectorFieldFunction CLinearDashpotInteraction::getVectorFieldFunction(const string& name)
{
CLinearDashpotInteraction::VectorFieldFunction vf;
if (name=="force"){
vf=&CLinearDashpotInteraction::getForce;
} else {
vf=NULL;
std::cerr << "ERROR - invalid name for interaction vector access function " << name << " in LinearDashpotInteraction" << endl;
}
return vf;
}
/*!
dummy
*/
CLinearDashpotInteraction::CheckedScalarFieldFunction CLinearDashpotInteraction::getCheckedScalarFieldFunction(const string& name)
{
CLinearDashpotInteraction::CheckedScalarFieldFunction csf;
csf=NULL;
cerr << "ERROR - invalid name for interaction vector access function " << name << " in LinearDashpotInteraction" << endl;
return csf;
}
|