← Back to branch summary

~ubuntu-branches/ubuntu/trusty/yade/trusty

« back to all changes in this revision

Viewing changes to pkg/dem/WirePM.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky, cf3f8d9
  • Date: 2013-10-30 20:56:33 UTC
  • mfrom: (20.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20131030205633-1f01r7hjce17d723
Tags: 1.05.0-2
http://bugs.debian.org/726009
[cf3f8d9] Pass -ftrack-macro-expansion=0 only if gcc>=4.8. (Closes: #726009)

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/dem/WirePM.cpp
10
10
#include<yade/core/Scene.hpp>
11
11
#include<yade/pkg/dem/ScGeom.hpp>
12
12
#include<yade/core/Omega.hpp>
 
13
#include "../../lib/base/Math.hpp"
 
14
 
 
15
#include <boost/random/linear_congruential.hpp>
 
16
#include <boost/random/triangle_distribution.hpp>
 
17
#include <boost/random/variate_generator.hpp>
 
18
 
 
19
#include<yade/core/Timing.hpp>
13
20
 
14
21
YADE_PLUGIN((WireMat)(WireState)(WirePhys)(Ip2_WireMat_WireMat_WirePhys)(Law2_ScGeom_WirePhys_WirePM));
15
22
 
22
29
        //BUG: ????? postLoad is called twice,
23
30
        LOG_TRACE( "WireMat::postLoad - update material parameters" );
24
31
        
25
 
        // compute cross-section area
 
32
        // compute cross-section area for single wire
26
33
        as = pow(diameter*0.5,2)*Mathr::PI;
27
34
 
 
35
        // check for stress strain curve for single wire
28
36
        if(strainStressValues.empty()) return; // uninitialized object, don't do nothing at all
29
37
        if(strainStressValues.size() < 2)
30
 
                throw invalid_argument("WireMat.strainStressValues: at least two points must be given.");
 
38
                throw std::invalid_argument("WireMat.strainStressValues: at least two points must be given.");
31
39
        if(strainStressValues[0](0) == 0. && strainStressValues[0](1) == 0.)
32
 
                throw invalid_argument("WireMat.strainStressValues: Definition must start with values greather then zero (strain>0,stress>0)");
 
40
                throw std::invalid_argument("WireMat.strainStressValues: Definition must start with values greater than zero (strain>0,stress>0)");
 
41
        
 
42
        switch(type) {
 
43
                case 0:
 
44
                        LOG_DEBUG("WireMat - Bertrand's approach");
 
45
                        if(!strainStressValuesDT.empty())
 
46
                                throw std::invalid_argument("Use of WireMat.strainStressValuesDT has no effect!");
 
47
                        break;
 
48
                case 1:
 
49
                        // check stress strain curve four double twist if type=1
 
50
                        LOG_DEBUG("WireMat - New approach with two curves");
 
51
                        if(isDoubleTwist) {
 
52
                                if(strainStressValuesDT.empty())
 
53
                                        throw runtime_error("WireMat.strainStressValuesDT not defined");
 
54
                                if(strainStressValuesDT.size() < 2)
 
55
                                        throw std::invalid_argument("WireMat.strainStressValuesDT: at least two points must be given.");
 
56
                                if(strainStressValuesDT[0](0) == 0. && strainStressValuesDT[0](1))
 
57
                                        throw std::invalid_argument("WireMat.strainStressValuesDT: Definition must start with values greater than zero (strain>0,stress>0)");
 
58
                        }
 
59
                        break;
 
60
                case 2:
 
61
                        // check stress strain curve four double twist if type=2
 
62
                        LOG_DEBUG("WireMat - New approach with two curves and initial shift");
 
63
                        if(isDoubleTwist) {
 
64
                                if(strainStressValuesDT.empty())
 
65
                                        throw runtime_error("WireMat.strainStressValuesDT not defined");
 
66
                                if(strainStressValuesDT.size() < 2)
 
67
                                        throw std::invalid_argument("WireMat.strainStressValuesDT: at least two points must be given.");
 
68
                                if(strainStressValuesDT[0](0) == 0. && strainStressValuesDT[0](1))
 
69
                                        throw std::invalid_argument("WireMat.strainStressValuesDT: Definition must start with values greater than zero (strain>0,stress>0)");
 
70
                        }
 
71
                        break;
 
72
                default:
 
73
                        throw std::invalid_argument("WireMat.type: Type must be 0, 1 or 2.");
 
74
                        break;
 
75
        }
33
76
 
34
77
}
35
78
 
53
96
        /* get reference to values since values are updated/changed in order to take unloading into account */
54
97
        vector<Vector2r> &DFValues = phys->displForceValues;
55
98
        vector<Real> &kValues = phys->stiffnessValues;
 
99
        Real kn = phys->kn;
56
100
 
57
101
        Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance
58
102
 
74
118
        
75
119
        /* compute normal force Fn */
76
120
        Real Fn = 0.;
77
 
        if ( D > DFValues[0](0) )
78
 
                Fn = kValues[0] * (D-phys->plastD);
79
 
        else {
80
 
                bool isDone = false;
81
 
                unsigned int i = 0;
82
 
                while (!isDone && i < DFValues.size()) { 
83
 
                        i++;
 
121
        if ( D > DFValues[0](0) ) { // unloading
 
122
                LOG_TRACE("WirePM: Unloading");
 
123
                Fn = kn * (D-phys->plastD);
 
124
        }
 
125
        else { // loading
 
126
                LOG_TRACE("WirePM: Loading");
 
127
                for (unsigned int i=1; i<DFValues.size(); i++) { 
84
128
                        if ( D > DFValues[i](0) ) {
85
 
                                Fn = DFValues[i-1](1) + (D-DFValues[i-1](0))*kValues[i];
86
 
                                phys->plastD = D - Fn/kValues[0];
87
 
                                isDone = true;
 
129
                                Fn = DFValues[i-1](1) + (D-DFValues[i-1](0))*kValues[i-1];
 
130
                                phys->plastD = D - Fn/kn;
88
131
                                // update values for unloading
89
132
                                DFValues[0](0) = D;
90
133
                                DFValues[0](1) = Fn;
 
134
                                break;
91
135
                        }
92
136
                }
93
137
        }
94
138
        
95
 
        TRVAR3( displN, D, Fn );
96
 
 
97
139
        /* compression forces cannot be applied to wires */
98
140
        if (Fn > 0.) Fn = 0.;
99
141
        
 
142
        TRVAR3( displN, D, Fn );
 
143
        
100
144
        phys->normalForce = Fn*geom->normal; // NOTE: normal is position2-position1 - It is directed from particle1 to particle2
101
145
 
102
146
        /* compute a limit value to check how far the interaction is from failing */
118
162
        }
119
163
        
120
164
        /* set shear force to zero */
121
 
        phys->shearForce = Vector3r::Zero();;
 
165
        phys->shearForce = Vector3r::Zero();
122
166
 
123
167
}
124
168
 
129
173
        
130
174
        /* avoid any updates if interactions which already exist */
131
175
        if(interaction->phys) return; 
 
176
        //TODO: make boolean to make sure physics are never updated, optimisation of contact detection mesh (no contact detection after link is created) 
132
177
        
133
178
        LOG_TRACE( "Ip2_WireMat_WireMat_WirePhys::go - create interaction physics" );
134
179
        
137
182
 
138
183
        /* set equilibrium distance, e.g. initial distance between particle (stress free state) */
139
184
        shared_ptr<WirePhys> contactPhysics(new WirePhys()); 
140
 
        contactPhysics->initD = geom->penetrationDepth;
 
185
        Real initD = geom->penetrationDepth;
141
186
        contactPhysics->normalForce = Vector3r::Zero();
142
187
 
143
188
        /* get values from material */
151
196
        if ( mat1->id == mat2->id ) { // interaction of two bodies of the same material
152
197
                crossSection = mat1->as;
153
198
                SSValues = mat1->strainStressValues;
154
 
                if ( (mat1->isDoubleTwist) && (abs(interaction->getId1()-interaction->getId2())==1) ) // bodies which id differs by 1 are double twisted
 
199
                if ( (mat1->isDoubleTwist) && (abs(interaction->getId1()-interaction->getId2())==1) ) {// bodies which id differs by 1 are double twisted
155
200
                        contactPhysics->isDoubleTwist = true;
156
 
                else
 
201
                        if ( mat1->type==1 || mat1->type==2 ) {
 
202
                                SSValues = mat1->strainStressValuesDT;
 
203
                                crossSection *= 2.;
 
204
                        }
 
205
                }
 
206
                else {
157
207
                        contactPhysics->isDoubleTwist = false;
 
208
                }
158
209
        }
159
210
        else { // interaction of two bodies of two different materials, take weaker material and no double-twist
160
211
                contactPhysics->isDoubleTwist = false;
168
219
                }
169
220
        }
170
221
        
171
 
        if(SSValues.empty()) throw std::invalid_argument("WireMat.strainStressValue is empty!");
172
 
        
173
222
        Real R1 = geom->radius1;
174
223
        Real R2 = geom->radius2;
175
224
        
176
 
        Real l0 = R1 + R2 - contactPhysics->initD; // initial length of the wire (can be single or double twisted)
 
225
        Real l0 = R1 + R2 - initD; // initial length of the wire (can be single or double twisted)
 
226
 
 
227
        /* compute displacement-force values */
 
228
        vector<Vector2r> DFValues;
 
229
        vector<Real> kValues;
 
230
        Real dl = 0.;
 
231
        bool isShifted = false;
 
232
        
 
233
        /* account for random distortion if type=2 */
 
234
        if ( mat1->type==2 ) {
 
235
                isShifted = true;
 
236
                if (mat1->seed==-1)
 
237
                        dl = l0*mat1->lambdau;
 
238
                else {
 
239
                        // initialize random number generator
 
240
                        static boost::minstd_rand randGen(mat1->seed!=0?mat1->seed:(int)TimingInfo::getNow(true));
 
241
                        static boost::variate_generator<boost::minstd_rand&, boost::triangle_distribution<Real> > rnd(randGen, boost::triangle_distribution<Real>(0,0.5,1));
 
242
                        Real rndu = rnd();
 
243
                        TRVAR1( rndu );
 
244
                        dl = l0*mat1->lambdau*rndu;
 
245
                        isShifted = true;
 
246
                }
 
247
        }
 
248
        else if ( mat2->type==2 ) {
 
249
                isShifted = true;
 
250
                if (mat2->seed==-1)
 
251
                        dl = l0*mat2->lambdau;
 
252
                else {
 
253
                        // initialize random number generator
 
254
                        static boost::minstd_rand randGen(mat2->seed!=0?mat2->seed:(int)TimingInfo::getNow(true));
 
255
                        static boost::variate_generator<boost::minstd_rand&, boost::triangle_distribution<Real> > rnd(randGen, boost::triangle_distribution<Real>(0,0.5,1));
 
256
                        Real rndu = rnd();
 
257
                        TRVAR1( rndu );
 
258
                        dl = l0*mat2->lambdau*rndu;
 
259
                }
 
260
        }
 
261
        contactPhysics->dL=dl;
 
262
        contactPhysics->isShifted=isShifted;
 
263
        
 
264
        // update geometry values
 
265
        l0 += dl;
 
266
        contactPhysics->initD = initD;
177
267
 
178
268
        /* compute threshold displacement-force values (tension negative since ScGem is used!) */
179
 
        vector<Vector2r> DFValues;
180
269
        for ( vector<Vector2r>::iterator it = SSValues.begin(); it != SSValues.end(); it++ ) {
181
270
                Vector2r values = Vector2r::Zero();
182
 
                values(0) = -(*it)(0)*l0;
 
271
//              values(0) = -(*it)(0)*l0;
 
272
                values(0) = -(*it)(0)*l0-dl;
183
273
                values(1) = -(*it)(1)*crossSection;
184
274
                DFValues.push_back(values);
185
275
        }
186
276
 
187
 
        /* compute elastic stiffness of single wire */
188
 
        vector<Real> kValues;
189
 
        Real k = DFValues[0](1) / DFValues[0](0);
 
277
        /* compute elastic stiffness for unloading*/
 
278
        Real k = DFValues[0](1) / (DFValues[0](0)+dl);
190
279
 
191
 
        /* update values if the interaction is double twisted */
192
 
        if ( contactPhysics->isDoubleTwist ) {
 
280
        /* update values if the interaction is a double twist and type=0 */
 
281
        if ( contactPhysics->isDoubleTwist && mat1->type==0 ) {
 
282
                // type=0 (force displacement values are computed by manipulating the values of the single wire by using the parameters lambdak and lambdaEps)
193
283
                Real alpha = atan( l0 / (3.*Mathr::PI*mat1->diameter) );
194
284
                Real kh = k * ( l0*mat1->diameter/crossSection ) / ( 48.*cos(alpha) * ( 41./9.*(1.+mat1->poisson) + 17./4.*pow(tan(alpha),2) ) );
195
285
                k = 2. * ( mat1->lambdak*kh + (1-mat1->lambdak)*k );
201
291
                        DFValues[i](1) *= mappingF;
202
292
                }
203
293
        }
 
294
        else {
 
295
        // type=1 and type=2 (force displacement values have already been computed by given stress-strain curve)
 
296
        } 
204
297
        
205
 
        /* store displacement-force values in physics */
206
 
        contactPhysics->displForceValues = DFValues;
 
298
        /* store elastic/unloading stiffness as kn in physics */
 
299
        contactPhysics->kn = k;
 
300
        contactPhysics->ks = 0.;
 
301
        TRVAR1( k );
 
302
 
 
303
        /* consider an additional point for the initial shift if type==2 */
 
304
        if ( mat1->type==2 ) {
 
305
                Vector2r values = Vector2r::Zero();
 
306
                values(0) = -dl+mat1->lambdaF*(DFValues[0](0)+dl);
 
307
                values(1) = DFValues[0](1)*mat1->lambdaF;
 
308
                k = values(1) / values(0);
 
309
                if ( mat1->lambdaF<1. )
 
310
                        DFValues.insert( DFValues.begin(), values );
 
311
        }
 
312
        else if ( mat2->type==2 ) {
 
313
                Vector2r values = Vector2r::Zero();
 
314
                values(0) = -dl+mat2->lambdaF*(DFValues[0](0)+dl);
 
315
                values(1) = DFValues[0](1)*mat2->lambdaF;
 
316
                k = values(1) / values(0);
 
317
                if ( mat2->lambdaF<1. )
 
318
                        DFValues.insert( DFValues.begin(), values );
 
319
        }
207
320
 
208
321
        /* compute stiffness-values of wire */
209
 
        contactPhysics->kn = k;
210
322
        kValues.push_back(k);
211
323
        for( unsigned int i = 1 ; i < DFValues.size(); i++ ) {
212
324
                Real deltau = -DFValues[i](0) + DFValues[i-1](0);
214
326
                k = deltaF/deltau;
215
327
                kValues.push_back(k);
216
328
        }
 
329
        
 
330
        /* add zero values for first point */
 
331
        DFValues.insert( DFValues.begin(), Vector2r::Zero() );
 
332
 
 
333
        /* store values in physics */
 
334
        contactPhysics->displForceValues = DFValues;
217
335
        contactPhysics->stiffnessValues = kValues;
218
336
 
219
337
        /* set particles as linked */
227
345
}
228
346
 
229
347
WirePhys::~WirePhys(){}
 
348