« back to all changes in this revision
Viewing changes to Model/ARotBondedInteraction.cpp
- browse files at revision 1204
- compare with another revision
- download diff
- download tarball
- view history from revision 1204
- Committer: Steffen Abe
- Date: 2019-07-30 11:28:05 UTC
- Revision ID: s.abe@igem-energie.de-20190730112805-73rdmsjtdqy17bnw
- EXPERIMENTAL added new brittle beam interaction type. N.B. untested,
should not be used for production!
- fixed cmake issue w.r.t. esys/lsm/geometry
should not be used for production!
- fixed cmake issue w.r.t. esys/lsm/geometry
- files added:
- Model/ARotBondedInteraction.cpp
- files modified:
- Foundation/Matrix3.cpp
added
removed
Model/ARotBondedInteraction.cpp
1
/////////////////////////////////////////////////////////////
2
// //
3
// Copyright (c) 2003-2017 by The University of Queensland //
4
// Centre for Geoscience Computing //
5
// http://earth.uq.edu.au/centre-geoscience-computing //
6
// //
7
// Primary Business: Brisbane, Queensland, Australia //
8
// Licensed under the Open Software License version 3.0 //
9
// http://www.apache.org/licenses/LICENSE-2.0 //
10
// //
11
/////////////////////////////////////////////////////////////
12
13
#include "Model/ARotBondedInteraction.h"
14
#include "Model/BondedInteractionCpData.h"
15
#include "Foundation/console.h"
16
17
// --- TML includes ---
18
#include "tml/message/packed_message_interface.h"
19
20
#include <stdexcept>
21
22
double calc_angle( double s_in, double c_os)
23
{
24
double angle;
25
if (s_in >0.0)
26
{
27
if ( c_os>=1.0 ||c_os<=-1.0)
28
{
29
angle = 0.0;
30
}
31
else {
32
angle = acos(c_os);
33
}
34
}
35
else if (s_in==0.0)
36
{
37
angle = 0.0;
38
}
39
else
40
{
41
if ( c_os>=1.0 ||c_os<=-1.0)
42
{
43
angle = 0.0;
44
}
45
else
46
{
47
angle = - acos(c_os);
48
}
49
}
50
return angle ;
51
}
52
53
54
/*!
55
Default constructor for basic rotational bonded interaction parameters.
56
Sets all stiffnesses to 0.0
57
*/
58
ARotBondedIGP::ARotBondedIGP()
59
: AIGParam(),
60
kr(0.0),
61
ks(0.0),
62
kt(0.0),
63
kb(0.0),
64
tag(0),
65
scaling(true)
66
{}
67
68
/*!
69
Constructor for basic rotational bonded interaction parameters from
70
individual stiffness parameters (Wang et al 2006)
71
*/
72
ARotBondedIGP::ARotBondedIGP(
73
const std::string &name,
74
double kr,
75
double ks,
76
double kt,
77
double kb,
78
int tag,
79
bool scaling
80
)
81
: AIGParam(name),
82
kr(kr),
83
ks(ks),
84
kt(kt),
85
kb(kb),
86
tag(tag),
87
scaling(scaling)
88
{}
89
90
/*!
91
Constructor for basic rotational bonded interaction parameters from
92
brittle beam parameters (stiffness only, not strength).
93
*/
94
ARotBondedIGP::ARotBondedIGP(
95
const std::string &name,
96
double youngsModulus,
97
double poissonsRatio,
98
int tag
99
)
100
: AIGParam(name),
101
tag(tag),
102
scaling(true)
103
{
104
double shearModulus = youngsModulus / (2.0*(1.+poissonsRatio));
105
106
kr = M_PI*youngsModulus;
107
ks = M_PI*shearModulus;
108
kb = M_PI*youngsModulus;
109
kt = M_PI*shearModulus;
110
}
111
112
// --------- END PARAMETER CLASS ----------
113
114
115
/*!
116
Default constructor for base of rotational bonded interactions.
117
Just sets all parameters to 0.0
118
*/
119
ARotBondedInteraction::ARotBondedInteraction():ARotPairInteraction()
120
{
121
m_kr = 0.0 ;
122
m_ks = 0.0 ;
123
m_kb = 0.0 ;
124
m_kt = 0.0 ;
125
126
m_nForce = 0.0 ;
127
m_shForce = 0.0;
128
m_tMoment = 0.0;
129
m_bMoment = 0.0;
130
m_moment = Vec3(0.0,0.0,0.0);
131
m_force = Vec3(0.0,0.0,0.0);
132
133
m_tag = 0;
134
m_scaling = true;
135
}
136
137
/*!
138
Constructor for base of rotational bonded interactions using two particle pointers
139
and a parameter set contained in a ARotBondedIGP.
140
*/
141
ARotBondedInteraction::ARotBondedInteraction(CRotParticle* p1,CRotParticle* p2,const ARotBondedIGP& param):ARotPairInteraction(p1,p2)
142
{
143
m_nForce = 0.0;
144
m_shForce = 0.0;
145
m_tMoment = 0.0;
146
m_bMoment = 0.0;
147
148
m_D = p2->getPos()-p1->getPos();
149
m_dist=sqrt(m_D*m_D);
150
m_r0=m_dist;
151
152
double effR=1.0;
153
double effL=1.0;
154
double effA=1.0;
155
double momentJ=1.0;
156
double momentI=1.0;
157
m_scaling = param.scaling;
158
159
// scale elastic param
160
// assume Rmean-scaling -> add option to do otherwise to derived class
161
if (m_scaling == true) {
162
if(!CParticle::getDo2dCalculations()){
163
effR=0.5*m_dist;
164
effL = m_dist;
165
effA = effR * effR;
166
momentJ = 0.5 * effR * effR * effR * effR;
167
momentI = 0.5 * momentJ;
168
}
169
}
170
171
m_kr = param.kr * effA / effL;
172
m_ks = param.ks * effA / effL;
173
m_kb = param.kb * momentI / effL;
174
m_kt = param.kt * momentJ / effL;
175
176
m_force = Vec3(0.0,0.0,0.0);
177
m_moment = Vec3(0.0,0.0,0.0) ;
178
179
m_tag = param.tag;
180
}
181
182
/*!
183
Default destructor - does nothing
184
*/
185
ARotBondedInteraction::~ARotBondedInteraction()
186
{}
187
188
/*!
189
get bond tag
190
*/
191
int ARotBondedInteraction::getTag() const
192
{
193
return m_tag;
194
}
195
196
/*!
197
set bond tag
198
*/
199
void ARotBondedInteraction::setTag(int tag)
200
{
201
m_tag = tag;
202
}
203
204
/*!
205
Get _initial_ vector between particle centers (i.e. at contruction time)
206
*/
207
Vec3 ARotBondedInteraction::getInitialCentrePtDiff() const
208
{
209
return m_p2->getInitPos() - m_p1->getInitPos();
210
}
211
212
/*!
213
Get _current_ vector between particle centers
214
*/
215
Vec3 ARotBondedInteraction::getCentrePtDiff() const
216
{
217
return m_p2->getPos() - m_p1->getPos();
218
}
219
220
/*!
221
Get location of the midpoint between the particles at construction time
222
*/
223
Vec3 ARotBondedInteraction::getInitialMidPoint() const
224
{
225
const Vec3 initialDiff = getInitialCentrePtDiff();
226
const Vec3 normalisedDiff = initialDiff*(1.0/initialDiff.norm());
227
return 0.5*((m_p1->getInitPos() + normalisedDiff*m_p1->getRad()) + (m_p2->getInitPos() - normalisedDiff*m_p2->getRad()));
228
}
229
230
/*!
231
Get shear force on particle 2
232
*/
233
Vec3 ARotBondedInteraction::getP2ShearForcePt() const
234
{
235
return m_p2->getPos() + (m_p2->getQuat().to_matrix()*(getInitialMidPoint()-m_p2->getInitPos()));
236
}
237
238
/*!
239
Get shear force on particle 1
240
*/
241
Vec3 ARotBondedInteraction::getP1ShearForcePt() const
242
{
243
return m_p1->getPos() + (m_p1->getQuat().to_matrix()*(getInitialMidPoint()-m_p1->getInitPos()));
244
}
245
246
/*!
247
?
248
*/
249
Vec3 ARotBondedInteraction::getShearDiff() const
250
{
251
const Vec3 p1Pt = getP1ShearForcePt();
252
const Vec3 p2Pt = getP2ShearForcePt();
253
const Vec3 ptDiff = p2Pt-p1Pt;
254
const Vec3 rotatedInitialDiffDirection = p1Pt - m_p1->getPos(); //getCentrePtDiff();
255
return (ptDiff - (((ptDiff*rotatedInitialDiffDirection)/rotatedInitialDiffDirection.norm2())*rotatedInitialDiffDirection));
256
}
257
258
/*!
259
get _current_ midpoint between particles
260
*/
261
Vec3 ARotBondedInteraction::getContactPoint() const
262
{
263
const Vec3 centrePtDiff = getCentrePtDiff();
264
const Vec3 normalisedDiff = centrePtDiff*(1.0/centrePtDiff.norm());
265
return 0.5*((m_p1->getPos() + normalisedDiff*m_p1->getRad()) + (m_p2->getPos() - normalisedDiff*m_p2->getRad()));
266
}
267
268
#if 0
269
270
Vec3 getTangent(const Vec3 &r, const Vec3 &f)
271
{
272
return (f.norm2() > 0) ? (f - ((r*f)/f.norm2())*r) : Vec3::ZERO;
273
}
274
275
void CRotBondedInteraction::calcForces()
276
{
277
// Calculate linear elastic force, no moment
278
const Vec3 currentDiff = m_p2->getPos() - m_p1->getPos();
279
const double currentDiffNorm = currentDiff.norm();
280
const Vec3 initialDiff = getInitialCentrePtDiff();
281
const Vec3 p2LinearForce = ((m_kr * (initialDiff.norm() - currentDiffNorm))/currentDiffNorm) * currentDiff;
282
283
m_p2->applyForce( p2LinearForce, m_p2->getPos());
284
m_p1->applyForce(-p2LinearForce, m_p1->getPos());
285
m_nForce = p2LinearForce.norm();
286
287
// Calculate the shearing forces, linear and moment contributions.
288
// Need to take into account relative orientations of particles.
289
const Vec3 shearDiff = getShearDiff();
290
const Vec3 p2ShearForcePt = getP2ShearForcePt();
291
const Vec3 p1ShearForcePt = getP1ShearForcePt();
292
const Vec3 p2ShearForce = -m_ks * (p2ShearForcePt - p1ShearForcePt); //(m_ks * shearDiff);
293
m_p2->applyMoment(cross(p2ShearForcePt - m_p2->getPos(), getTangent(p2ShearForcePt - m_p2->getPos(), p2ShearForce)));
294
m_p1->applyMoment(cross(p1ShearForcePt - m_p1->getPos(), getTangent(p1ShearForcePt - m_p1->getPos(), -p2ShearForce)));
295
m_p2->applyForce( p2ShearForce, m_p2->getPos());
296
m_p1->applyForce(-p2ShearForce, m_p1->getPos());
297
m_shForce = p2ShearForce.norm();
298
}
299
300
#else
301
302
const double HALF_SQRT_2 = 0.5*sqrt(2.0); // use constexpr?
303
304
/*!
305
Calculate forces
306
*/
307
void ARotBondedInteraction::calcForces()
308
{
309
double s_fai=0.0, c_fai=0.0, c_pasi2=0.0, s_pasi2=0.0;
310
double sita=0.0, pasi=0.0;
311
312
const Matrix3 mat0 = (m_p1->getQuat()).to_matrix();
313
const Vec3 rbp = m_p2->getPos() - m_p1->getPos(); // rbp <-> r_f
314
const double rbpNorm = rbp.norm();
315
const Vec3 rb = mat0*(m_p2->getPos() - m_p1->getPos()); // eq. 11, rb <-> r_c
316
const double rbNorm = rb.norm();
317
318
const double r_0Norm = m_p1->getRad()+m_p2->getRad();
319
//const double r_0Norm = m_r0;
320
const Vec3 delta_r = (rbNorm - r_0Norm)*rb/rbNorm; // eq. 9, part of eq. 12
321
const Vec3 Fr = m_kr*delta_r ; // rest of eq. 12
322
323
const double gama_cos = dot(rb, m_D)/(rbNorm*m_D.norm()) ;
324
double gama = 0.0;
325
if (gama_cos < 1.0 && gama_cos > -1.0 )
326
{
327
gama = acos(gama_cos);
328
}
329
330
const Vec3 rb_0 = cross(rb, m_D);
331
const Vec3 rb_b0 = cross(rb, rb_0);
332
const double rb_b0Norm = rb_b0.norm();
333
Vec3 Fst(0.0,0.0,0.0);
334
if (rb_b0Norm != 0)
335
{
336
Fst = m_ks*r_0Norm*gama*rb_b0/rb_b0Norm; // eq. 14
337
}
338
339
Vec3 Mst(0.0,0.0,0.0);
340
const double rb_0Norm = rb_0.norm();
341
if (rb_0Norm != 0)
342
{
343
Mst = -Fst.norm()*rb_0/rb_0Norm; // eq. 15
344
}
345
346
const double m0 = HALF_SQRT_2*sqrt((rbNorm + rb.Z())/rbNorm);
347
double m1 = -HALF_SQRT_2*sqrt((rbNorm - rb.Z())/rbNorm) ;
348
double m2 = -m1;
349
const double m3 = 0.0;
350
const double rbX2PlusRbY2 = rb.X()*rb.X() + rb.Y()*rb.Y();
351
if (rbX2PlusRbY2 != 0.0)
352
{
353
const double denomSqrt = sqrt(rbX2PlusRbY2);
354
m1 = m1*rb.Y()/denomSqrt;
355
m2 = m2*rb.X()/denomSqrt;
356
}
357
358
const Quaternion qm(m0, Vec3(m1,m2,m3));
359
const Matrix3 mat = qm.to_matrix();
360
const Matrix3 matTrans(mat.trans());
361
const Quaternion rp = (m_p1->getQuat()).inverse() * (m_p2->getQuat());
362
const Quaternion r = qm.inverse() * rp * qm;
363
const double temp0 = r.return_sca()*r.return_sca() +
364
r.return_vec().Z()*r.return_vec().Z();
365
if(temp0 == 0.0)
366
{
367
pasi = 0.0;
368
}
369
else
370
{
371
const double sqrtTemp0 = sqrt(temp0);
372
c_pasi2 = r.return_sca()/sqrtTemp0;
373
s_pasi2 = r.return_vec().Z()/sqrtTemp0;
374
pasi = 2.0 * calc_angle(s_pasi2,c_pasi2);
375
}
376
const Vec3 Mtp = Vec3(0.0,0.0, m_kt*pasi);
377
const Vec3 Mtr = matTrans * Mtp;
378
const double temp2 = r.return_vec().X()*r.return_vec().X() +
379
r.return_vec().Y()*r.return_vec().Y();
380
const double temp3 = r.return_vec().X()*r.return_vec().Z() +
381
r.return_sca() * r.return_vec().Y();
382
const double temp4 = r.return_vec().Y()*r.return_vec().Z() -
383
r.return_sca() * r.return_vec().X();
384
385
Vec3 Mbr(0, 0, 0);
386
Vec3 Fsr(0, 0, 0);
387
Vec3 Msr(0, 0, 0);
388
if (temp2 != 0.0)
389
{
390
const double temp1 = temp0 - temp2;
391
if (temp1 >= 1.0 || temp1 <= -1.0)
392
{
393
sita = 0.0 ;
394
}
395
else
396
{
397
sita = acos(temp1);
398
}
399
if (temp0 == 0.0) {
400
const double sqrtDenom = sqrt(temp2);
401
c_fai = r.return_vec().Y()/sqrtDenom;
402
s_fai = -r.return_vec().X()/sqrtDenom;
403
} else {
404
const double sqrtDenom = sqrt(temp0*temp2);
405
c_fai = temp3/sqrtDenom;
406
s_fai = temp4/sqrtDenom;
407
}
408
const Vec3 Mbp = Vec3( - 0.5*m_kb*sita*s_fai, 0.5*m_kb*sita*c_fai, 0.0); // eq 16.1
409
Mbr = matTrans * Mbp; // eq. 17.1
410
const Vec3 Fsp = Vec3 ( - 0.5*m_ks*c_fai*sita*r_0Norm,
411
- 0.5*m_ks*s_fai*sita*r_0Norm,
412
0.0); // eq. 16.3
413
Fsr = matTrans * Fsp; // eq. 17.3
414
415
const Vec3 Msp = Vec3( 0.5*m_ks*s_fai*sita*r_0Norm,
416
- 0.5*m_ks*c_fai*sita*r_0Norm,
417
0.0 );
418
// above is diff to paper: r instead r^2 here, other r moved to eq 19 below
419
Msr = matTrans * Msp;
420
}
421
422
const double eq_rad1 = m_p1->getRad()/(m_p1->getRad()+m_p2->getRad());
423
const double eq_rad2 = m_p2->getRad()/(m_p1->getRad()+m_p2->getRad());
424
425
const Matrix3 mat0Trans(mat0.trans());
426
m_force = mat0Trans * (Fr + Fst + Fsr );
427
m_moment = mat0Trans * (Mbr + Mtr + Msr*eq_rad1*rbpNorm + Mst*eq_rad1*rbpNorm); // eq. 19 for particle 1
428
429
const Vec3 moment2 = mat0Trans *( Msr*eq_rad2*rbpNorm + Mst*eq_rad2*rbpNorm - Mbr - Mtr); // eq. 19 for particle 2
430
431
m_shForce = (Fst + Fsr).norm();
432
m_tMoment = Mtr.norm();
433
m_bMoment = Mbr.norm();
434
const double r0 = m_p1->getRad()+m_p2->getRad();
435
const Vec3 D = m_p1->getPos() - m_p2->getPos();
436
m_dist = sqrt(D*D);
437
438
m_nForce = ((m_dist < r0) ? -1.0 : 1.0) * Fr.norm();
439
440
Vec3 pos=m_p2->getPos()+(m_p2->getRad()/(m_p1->getRad()+m_p2->getRad()))*D;
441
442
m_p2->applyForce(-1.0*m_force,pos);
443
m_p1->applyForce(m_force,pos);
444
m_cpos=pos;
445
446
m_p1->applyMoment(m_moment) ;
447
m_p2->applyMoment(moment2) ;
448
}
449
450
#endif
451
452
// --- FIELD FUNCTIONS ---
453
/*!
454
get stored elastic energy
455
*/
456
double ARotBondedInteraction::getPotentialEnergy() const
457
{
458
double pe_r,pe_s,pe_t,pe_b;
459
460
pe_r=(m_kr!=0.0) ? 0.5*(m_nForce*m_nForce)/m_kr : 0.0;
461
pe_s=(m_ks!=0.0) ? 0.5*(m_shForce*m_shForce)/m_ks : 0.0;
462
pe_t=(m_kt!=0.0) ? 0.5*(m_tMoment*m_tMoment)/m_kt : 0.0;
463
pe_b=(m_kb!=0.0) ? 0.5*(m_bMoment*m_bMoment)/m_kb : 0.0;
464
return pe_r+pe_s+pe_t+pe_b;
465
}
466
467
double ARotBondedInteraction::getNormalPotentialEnergy() const
468
{
469
return (m_kr!=0.0) ? 0.5*(m_nForce*m_nForce)/m_kr : 0.0;
470
}
471
472
double ARotBondedInteraction::getShearPotentialEnergy() const
473
{
474
return (m_ks!=0.0) ? 0.5*(m_shForce*m_shForce)/m_ks : 0.0;
475
}
476
477
double ARotBondedInteraction::getTwistPotentialEnergy() const
478
{
479
return (m_kt!=0.0) ? 0.5*(m_tMoment*m_tMoment)/m_kt : 0.0;
480
}
481
482
double ARotBondedInteraction::getBendPotentialEnergy() const
483
{
484
return (m_kb!=0.0) ? 0.5*(m_bMoment*m_bMoment)/m_kb : 0.0;
485
}
486
487
Vec3 ARotBondedInteraction::getForce() const
488
{
489
return m_force;
490
}
491
492
Vec3 ARotBondedInteraction::getNormalForce() const
493
{
494
Vec3 dr = (m_p1->getPos() - m_p2->getPos()).unit();
495
Vec3 normalForce = (m_force*dr)*dr;
496
return normalForce;
497
}
498
499
Vec3 ARotBondedInteraction::getTangentialForce() const
500
{
501
return (m_force - getNormalForce());
502
}
503
504
/*!
505
Get the particle member function which returns a scalar field of a given name.
506
507
\param name the name of the field
508
*/
509
ARotBondedInteraction::ScalarFieldFunction ARotBondedInteraction::getScalarFieldFunction(const string& name)
510
{
511
ARotBondedInteraction::ScalarFieldFunction sf;
512
513
if (name=="potential_energy"){
514
sf=&ARotBondedInteraction::getPotentialEnergy;
515
} else if (name=="e_pot_normal"){
516
sf=&ARotBondedInteraction::getNormalPotentialEnergy;
517
} else if (name=="e_pot_shear"){
518
sf=&ARotBondedInteraction::getShearPotentialEnergy;
519
} else if (name=="e_pot_twist"){
520
sf=&ARotBondedInteraction::getTwistPotentialEnergy;
521
} else if (name=="e_pot_bend"){
522
sf=&ARotBondedInteraction::getBendPotentialEnergy;
523
} else if (name=="count"){
524
sf=&ARotBondedInteraction::Count;
525
} else {
526
sf=NULL;
527
cerr << "ERROR - invalid name for interaction scalar access function" << endl;
528
}
529
530
return sf;
531
}
532
533
/*!
534
Get the particle member function which returns a vector field of a given name.
535
536
\param name the name of the field
537
*/
538
ARotBondedInteraction::VectorFieldFunction ARotBondedInteraction::getVectorFieldFunction(const string& name)
539
{
540
ARotBondedInteraction::VectorFieldFunction vf;
541
542
if (name=="force"){
543
vf=&ARotBondedInteraction::getForce;
544
} else if (name=="normal_force"){
545
vf=&ARotBondedInteraction::getNormalForce;
546
} else if (name=="tangential_force"){
547
vf=&ARotBondedInteraction::getTangentialForce;
548
} else {
549
vf=NULL;
550
cerr << "ERROR - invalid name for interaction vector access function" << endl;
551
}
552
553
return vf;
554
}
555
556
/*!
557
Get the particle member function which returns a checked scalar field of a given name.
558
559
\param name the name of the field
560
*/
561
ARotBondedInteraction::CheckedScalarFieldFunction ARotBondedInteraction::getCheckedScalarFieldFunction(const string& name)
562
{
563
ARotBondedInteraction::CheckedScalarFieldFunction sf;
564
565
sf=NULL;
566
cerr << "ERROR - invalid name for interaction scalar access function" << endl;
567
568
return sf;
569
}
570
Loggerhead is a web-based interface for Breezy
Version: 2.0.1