File: beads_column.rst

package info (click to toggle)
siconos 4.3.1%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 82,496 kB
  • sloc: cpp: 159,693; ansic: 108,665; fortran: 33,248; python: 20,709; xml: 1,244; sh: 385; makefile: 226
file content (498 lines) | stat: -rw-r--r-- 20,964 bytes parent folder | download | duplicates (3)
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
.. _beads_column_example:

Tutorial : a column of three beads
==========================================

.. highlight:: c++

This tutorial deals with a column of beads, subjected to the
gravity. It introduces Lagrangian Linear systems with Lagrangian
Linear Relations and Newton impact law.

Building the Non-Smooth Dynamical System
----------------------------------------

As described in the figure below, we consider a ball of mass :math:`m` and radius :math:`R`, described by 3 generalized coordinates :math:`q=(z,x,\theta)`
The ball is subjected to the gravity :math:`g`. The system is also constituted by a rigid plane, defined by its position :math:`h` with respect
to the axis Oz. We assume that the position of the plane is fixed and not modeled as a dynamical system.

.. figure:: /figures/mechanics/BouncingBall/BouncingBall.*
   :width: 8cm
   :align: center

   fig 1: Coordinate system

The equation of motion of a ball is given by

.. math::

   M\ddot q = F_{ext}(t) + p \ \ with \ \ M=\left[\begin{array}{ccc}
   m &0 &0 \\
   0 & m & 0 \\
   0 & 0 & I
   \end{array}\right] \ \ , \ \ I = \frac{3}{5}mR^2 \ \ and \ \ F_{ext} = \left[\begin{array}{c}
   -m g \\
   0  \\
   0
   \end{array}\right]

with

* :math:`M` the inertia term, a :math:`n\times{}n` matrix.
* :math:`p` the force due to the non-smooth law, ie the reaction at impact. 
* :math:`F_{ext}(t):  \mathcal R \rightarrow \mathcal R^{n}` the given external force.

That fits with Lagrangian, Linear and Time-Invariant Dynamical System, represented by LagrangianLinearTIDS class (see \ref dsInSiconos).

Next, we will suppose that we have a column of "dsNumber" balls like the one above (mass and radius may be different, and if necessary, variables will be indexed by :math:`i`, the number of the ball). Each ball is governed by a linear system like the one written for a single ball and may be in contact with the balls above and below it.

Let us now start the writing of the input file. Like for the first tutorial, we create a new directory, multiBeads, and save the template given \ref tutGCtemplate "here" as multiBeads.cpp. 

We start by setting some parameters, like the number of balls, their initial positions and velocities and so on::

  // User-defined main parameters 
  unsigned int dsNumber = 10;      // the number of dynamical systems 
  unsigned int nDof = 3;           // degrees of freedom for beads
  double increment_position = 1;   // initial position increment from one DS to the following
  double increment_velocity = 0;   // initial velocity increment from one DS to the following
  double t0 = 0;                   // initial computation time
  double T = 10;                   // final computation time 
  double h = 0.005;                // time step
  double position_init = 10.5;     // initial position for lowest bead.
  double velocity_init = 0.0;      // initial velocity for lowest bead.
  double R = 0.1;                  // balls radius

Then, we define some initial conditions for the balls, and create the
corresponding Dynamical Systems, all of type Lagrangian, Linear and
Time Invariant.

All the systems are inserted in a container, a DynamicalSystemsSet,
named allDS. 

From now on, to simplify writing, we suppose that all
balls have the same mass, :math:`m = 1`, and the same radius, :math:`R=0.1`::

  // -------------------------
  // --- Dynamical systems --- 
  // -------------------------
  
  // mass matrix, set to identity
  SP::SiconosMatrix Mass = new SimpleMatrix(nDof,nDof);
  Mass->eye();
  (*Mass)(2,2) = 3.0/5*R*R;
  
  // -- Initial positions and velocities --
  // q0[i] and v0[i] correspond to position and velocity of ball i. 
  vector<SimpleVector *> q0;
  vector<SimpleVector *> v0;
  q0.resize(dsNumber,NULL);
  v0.resize(dsNumber,NULL);
  
  for (unsigned i = 0; i < dsNumber; ++i)
  {
      // Memory allocation for q0[i] and v0[i]
      q0[i] = new SimpleVector(nDof);
      v0[i] = new SimpleVector(nDof); 
      // set values
      (*(q0[i]))(0) = position_init;
      (*(v0[i]))(0) = velocity_init;
      // Create a new Lagrangian Linear Dynamical System, with q0] and v0[i] as initial conditions,
      // Mass as mass matrix and i as number of identification. 
      // The system is then inserted in allDS. 
      allDS.insert( new LagrangianLinearTIDS(i,nDof,*(q0[i]),*(v0[i]),*Mass));
      // Increment values for next system
      position_init+= increment_position;
      velocity_init+= increment_velocity;
  }
  
Next, it is necessary to define the external forces, the gravity, applied on each ball. According to :ref:`dsPlugins`, a plug-in function is available for those forces. (For details on plug-in functions, see :ref:`siconos_plugins`). Its signature (the type of its arguments) is given in DefaultPlugin.cpp. So we copy it in a new file, say BeadsPlugin.cpp, and we define an extern function, gravity.::
  
  const double m = 1; // bead mass
  const double g = 9.81; // gravity
  extern "C" void gravity(unsigned int sizeOfq, double time,	double * fExt, double *param)
  {
      // set fExt components to 0
      for (unsigned int i = 0; i < sizeOfq; i++)
      fExt[i] = 0.0;
      // apply gravity
      fExt[0] = -m*g;
  } 

Warning

* gravity must be an *extern "C"* function, and code is C, not C++. 
* the name of the plugin file, BeadsPlugin.cpp here, must be xxxPlugin.cpp, xxx being whatever you want. 


Now we have to say "use gravity from BeadsPlugin.cpp to compute the external forces of my systems."
This is done thanks to "setComputeFExtFunction" function, in multiBeads.cpp::

   // 	
   CheckInsertDS checkDS; 
   for (i=0;i<dsNumber;i++)
      {
        // Memory allocation for q0[i] and v0[i]
	q0[i] = new SimpleVector(nDof);
	v0[i] = new SimpleVector(nDof); 
        // set values
	(*(q0[i]))(0) = position_init;
	(*(v0[i]))(0) = velocity_init;
        // Create and insert in allDS a new Lagrangian Linear Dynamical System ...
	checkDS = allDS.insert(new LagrangianLinearTIDS(i,nDof,*(q0[i]),*(v0[i]),*Mass));
        // Note that we now use a CheckInsertDS object: checkDS.first is 
	// an iterator that points to the DS inserted above.
        // 
	// Set the external forces for the last created system. 
	(static_cast<LagrangianDS*>(*(checkDS.first)))->setComputeFExtFunction("BeadsPlugin.so", "gravity");
	// A cast is required, since allDS handles DynamicalSystem*, 
	// not LagrangianLinearTIDS*. 
	// Increment values for next system
	position_init+= increment_position;
	velocity_init+= increment_velocity;
      }

From this point, any call to the external forces of a system in allDS will result in a call to the function gravity defined in BeadsPlugin.cpp.

*Remark:* :math:`m` and :math:`R` are set inside the BeadsPlugin file
but it would also be possible, and maybe better, to pass them as
parameters in gravity function.
See \ref doc_usingPlugin for details on that option.

Ok, now DynamicalSystems are clearly defined and all saved in allDS. Let's turn our attention to Interactions. In the same way, they will be handled by a container, an InteractionsSet, named allInteractions. The potential interactions are the contacts between beads and the impact on the ground. Thus, for dsNumbers systems, there are dsNumbers-1 "bead-bead" Interactions plus one between the "bottom bead" and the floor.

We start with bead-floor Interaction: the ball at the bottom bounces on the rigid plane, introducing a constraint on the position of the ball, given by:
:math:`z-R-h\geq 0`.
To define an Interaction, it is first necessary to set some relations between local variables at contact and the global coordinates. 
Thus, as a local variables of the Interaction, we introduce :math:`y` as the distance between the ball and the floor and :math:`\lambda` as the multiplier that corresponds to 
the reaction at contact. Then the relation is written,

.. math::

   y =& Hq + b = [1 \ 0 \ 0] q - R - h \\
   p =& H^t \lambda

(next, we set h=0).

.. compound::

   Finally we need to define a non-smooth law to define the behavior of the ball at impact. 
   The unilateral constraint is such that

   .. math:: 0 \leq y \perp \lambda \geq 0

   completed with a Newton Impact law, for which we set the restitutive coefficient :math:`e` to 0.9: 

   .. math:: \textrm{if} \ y=0, \ \dot y(t^+) = -e \dot y(t^-)

   with :math:`t^+` and :math:`t^-` being post and pre-impact times.

The first Interaction can then be constructed::

  // -------------------
  // --- Interactions---
  // -------------------
  InteractionsSet allInteractions;
  // The total number of Interactions
  int interactionNumber = dsNumber;
  // Interaction first bead and floor
  // A set for the systems handles by the "current" Interaction
  DynamicalSystemsSet dsConcerned;
  // Only the "bottom" bead is concerned by this first Interaction,
  // therefore DynamicalSystem number 0.
  dsConcerned.insert(allDS.getDynamicalSystemPtr(0)); 
  // -- Newton impact law -- 
  double e = 0.9;
  NonSmoothLaw * nslaw0 = new NewtonImpactNSL(e);
  // Lagrangian Relation
  unsigned int interactionSize = 1; // y vector size
  SiconosMatrix *H = new SimpleMatrix(interactionSize,nDof);
  (*H)(0,0) = 1.0;
  SiconosVector *b = new SimpleVector(interactionSize);
  (*b)(0) = -R;
  Relation * relation0 = new LagrangianLinearR(*H,*b);
  // Interaction
  unsigned int num = 0 ; // an id number for the Interaction
  Interaction * inter0 = new Interaction("bead-floor", dsConcerned,num,interactionSize, nslaw0, relation0);
  allInteractions.insert(inter0);

In the same way, the potential contact between two balls introduces some new constraints:

:math:`(z_i-R_i)-(z_j-R_j)-h \geq 0`, if ball :math:`i` is on top of ball :math:`j`.

So if we consider the Interaction between ball :math:`i` and :math:`j`, :math:`y` being the distance between two balls and :math:`\lambda` the multiplier, we get:

.. math::

   y =& HQ + b = [-1 \ 0 \ 0 \ 1 \ 0 \ 0]Q + R_j-R_i-h \\
   p =& H^t \lambda \\
   Q =& \left[\begin{array}{c}
   q_j \\
   q_i
   \end{array}\right]

With the same non smooth law as for the first Interaction::

  // A list of names for the Interactions
  vector<string> id;
  id.resize(interactionNumber-1);
  CheckInsertInteraction checkInter;
  // A vector that will handle all the relations
  vector<Relation*> LLR(interactionNumber-1);
  // 
  SiconosMatrix *H1 = new SimpleMatrix(1,2*nDof);
  if (dsNumber>1)
  {
      (*H1)(0,0) = -1.0;
      (*H1)(0,3) = 1.0;
      // Since Ri=Rj and h=0, we do not need to set b.
      Relation * relation = new LagrangianLinearR(*H1);
      for (i=1;(int)i<interactionNumber;i++)
      {
          // The systems handled by the current Interaction ...
	  dsConcerned.clear();
	  dsConcerned.insert(allDS.getDynamicalSystemPtr(i-1));
	  dsConcerned.insert(allDS.getDynamicalSystemPtr(i));
	  // The id: "i"
	  ostringstream ostr;
	  ostr << i;
	  id[i-1]= ostr.str();
	  // The relations
	  LLR[i-1] = new LagrangianLinearR(*relation); // we use copy constructor to built all relations
	  checkInter = allInteractions.insert( new Interaction(id[i-1], dsConcerned,i,interactionSize, nslaw0, LLR[i-1]));
      }
      delete relation;
  }

Note that each Relation corresponds to one and only one Interaction (which is not the case of NonSmoothLaw); that's why we need to built a new Relation LLR[i-1] for each Interaction. 

Everything is now ready to build the NonSmoothDynamicalSystem and the related Model::

    // --------------------------------
    // --- NonSmoothDynamicalSystem --- 
    // --------------------------------
    NonSmoothDynamicalSystem * nsds = new NonSmoothDynamicalSystem(allDS, allInteractions);    
    // -------------
    // --- Model ---
    // -------------
    Model * multiBeads = new Model(t0,T); 
    multiBeads->setNonSmoothDynamicalSystemPtr(nsds); // set NonSmoothDynamicalSystem of this model

The Simulation
--------------

Time-Stepping scheme
""""""""""""""""""""

As a first example, we will use a Moreau's time-stepping scheme, where the non-smooth problem will be written as a LCP. The process is more or less the same as for the Diode Bridge case, so we won't detail it. The only difference is that now, the OneStepIntegrator handles several DynamicalSystems::

  string solverName = "Lemke";      // solver algorithm used for non-smooth problem
  Simulation* s = new TimeStepping(multiBeads);
  // -- Time discretisation --
  TimeDiscretisation * t = new TimeDiscretisation(h,s);
  // -- OneStepIntegrators --
  double theta = 0.5000001; 
  OneStepIntegrator * OSI = new Moreau(allDS , theta ,s);
  // That means that all systems in allDS have the same theta value.	 
  // -- OneStepNsProblem --
  OneStepNSProblem * osnspb = new LCP(s,"LCP",solverName,10001, 0.001);

Event-Driven algorithm
""""""""""""""""""""""

In that second part, an event-driven algorithm is used to solve the problem. Event-Driven Simulation principle is detailed in :ref:`event_driven`.

The dynamics is decomposed in "modes", time-intervalls where the dynamics is smooth and discrete events where the dynamics is non-smooth.

In the present case, non smooth events will corresponds to impacts between balls. Each time such an event is detected, a non-smooth problem is formalized and solved (as a LCP here) while between events, the systems are integrated thanks to Lsodar, ODE solver with roots-finding algorithm.

As for the Time-stepping, we first need to built the simulation and then its time-discretisation::

  // The simulation belongs to Model multiBeads
  EventDriven* s = new EventDriven(multiBeads); 
  TimeDiscretisation * t = new TimeDiscretisation(h,s);

Next step is the declaration of integrators for the dynamical systems.
The integrator will handle all the DynamicalSystems of the Model. During integration of the systems, Lsodar will search for roots of some equations (the constraints ie the Interactions of the NonSmoothDynamicalSystem). The required OSI type is Lsodar, applied to allDS::
  
  OneStepIntegrator * OSI = new Lsodar(allDS,s); 

Each time a root is found, a new NonSmoothEvent is created and it's then necessary to write and solve a non-smooth problem. We won't detail this here but just remember that this requires two LCP, one at "velocity" level, named impact, and another at "acceleration" level, named acceleration. 
The whole event-driven algorithm for Lagrangian Systems is available here: :ref:`event_driven_lagrange`::

  OneStepNSProblem * impact = new LCP(s, "impact",solverName,101, 0.0001,"max",0.6);
  OneStepNSProblem * acceleration = new LCP(s, "acceleration",solverName,101, 0.0001,"max",0.6);

The Model is now complete, we can start the simulation process.

Simulation Process
------------------

Time-Stepping
"""""""""""""

Once again, the process is the same as in the first tutorial and won't be detailed.
Concerning the output, we save the position and velocity of all balls::

  s->initialize(); 
  int k = 0;
  int N = t->getNSteps(); // Number of time steps
  // Prepare output and save value for the initial time
  unsigned int outputSize = dsNumber*2+1;
  SimpleMatrix dataPlot(N+1,outputSize ); // Output data matrix
  // time
  dataPlot(k, 0) = multiBeads->getT0();
  // Positions and velocities
  i = 0; // Remember that DS are sorted in a growing order according to their number.
  DSIterator it;
  for(it = allDS.begin();it!=allDS.end();++it)
  {
      dataPlot(k,(int)i*2+1) = static_cast<LagrangianLinearTIDS*>(*it)->getQ()(0);
      dataPlot(k,(int)i*2+2) = static_cast<LagrangianLinearTIDS*>(*it)->getVelocity()(0);
      i++;
  }

Note that we use a "DSIterator", which is simply a pointer to a set of DynamicalSystems; allDS.begin() is a pointer to the first object handled by allDS and allDS.end() a pointer "just after" the last object handled by allDS. The current pointed system is then \*it ("content of the pointer"). Thus, in the loop above, we sweep through all the DynamicalSystems and get the corresponding :math:`q` and :math:`v`.
A static_cast is also required since allDS contains DynamicalSystem whereas we need functions specific to LagrangianDS (getQ ...). 

Next, we write::

  while(k < N)	
  {
      k++;	
      // solve ... 
      s->computeOneStep();
      dataPlot(k, 0) = s->getNextTime();
      // 
      i = 0;
      for(it = allDS.begin();it!=allDS.end();++it)
      {
          dataPlot(k,(int)i*2+1) = static_cast<LagrangianLinearTIDS*>(*it)->getQ()(0);
	  dataPlot(k,(int)i*2+2) = static_cast<LagrangianLinearTIDS*>(*it)->getVelocity()(0);
	  i++;
	  s->nextStep();
      }
  }

and for output file saving::

  ioMatrix io("result.dat", "ascii");
  io.write(dataPlot,"noDim");

Event-Driven
""""""""""""

The principle of an EventDriven simulation roughly consists in integration between some events with stops and special treatment at these events. Thus we introduce a specific object, the EventsManager, a kind of stack of events used to handle them, where they are saved in a chronological order. It belongs to the Simulation object and can be accessed with::

  EventsManager * eventsManager = s->getEventsManagerPtr();

The manager is built during the ininitialization, which is still the first required step of any simulation process::

  s->initialize();

Among other things, this initialization schedules time events from the TimeDiscretisation object into the manager. Each time step is saved as a TimeDiscretionEvent.

Then the simulation process consists in:
* check if there is a "future" event
* integrate the system until this future event is reached or until a non-smooth event is found
* schedule the possibly new event
* deal with the system at event (for example, in case of a non-smooth event, formalize and solve one or more LCP)
* next step

Once again this is only a summary and we encourage you to read :ref:`event_driven` to get more details about the event-driven strategy. 

The resulting code is::

   // While there are some events in the manager ...
    while(eventsManager->hasNextEvent())
      {
	eventDriven->computeOneStep();
      }

Concerning output, we first save displacements and velocities at each time step::

    while(eventsManager->hasNextEvent())
      {
	k++;
	eventDriven->advanceToEvent();

	eventDriven->processEvents();
        // Positions and velocities for user time steps
	i = 0; // Remember that DS are sorted in a growing order according to their number.
	DSIterator it;
	dataPlot(k, 0) = eventDriven->getStartingTime(); 
	for(it = allDS.begin();it!=allDS.end();++it)
	  {
	    dataPlot(k,(int)i*2+1) = static_cast<LagrangianLinearTIDS*>(*it)->getQ()(0);
	    dataPlot(k,(int)i*2+2) = static_cast<LagrangianLinearTIDS*>(*it)->getVelocity()(0);
	    i++;
	  }
      }

But when a non-smooth event occurs, that may be interesting to get pre and post impact values. 
In Siconos, the values saved in object are usually the last computed, thus in the present case, post-impact values.
The next-to-last values are saved in "memory" objects; we get them in case of "Non-Smooth event"::

    while(eventsManager->hasNextEvent())
      {
	k++;
	eventDriven->advanceToEvent();

	eventDriven->processEvents();
	if(eventsManager->getStartingEventPtr()->getType() == "NonSmoothEvent")
	  {
	    i = 0; // Remember that DS are sorted in a growing order according to their number.
	    DSIterator it;
	    dataPlot(k, 0) = eventDriven->getStartingTime(); 
	    for(it = allDS.begin();it!=allDS.end();++it)
	      {
		dataPlot(k,(int)i*2+1) = (*static_cast<LagrangianLinearTIDS*>(*it)->getQMemoryPtr()->getSiconosVector(1))(0);
		dataPlot(k,(int)i*2+2) = (*static_cast<LagrangianLinearTIDS*>(*it)->getVelocityMemoryPtr()->getSiconosVector(1))(0);
		i++;
	      }
	    k++;
	  }
        // Positions and velocities for user time steps
	i = 0; // Remember that DS are sorted in a growing order according to their number.
	DSIterator it;
	dataPlot(k, 0) = eventDriven->getStartingTime();  
	for(it = allDS.begin();it!=allDS.end();++it)
	  {
	    dataPlot(k,(int)i*2+1) = static_cast<LagrangianLinearTIDS*>(*it)->getQ()(0);
	    dataPlot(k,(int)i*2+2) = static_cast<LagrangianLinearTIDS*>(*it)->getVelocity()(0);
	    i++;
	  }
      }

    // Output written in result.dat 
    ioMatrix io("result.dat", "ascii");
    io.write(dataPlot,"noDim");

The simulation is now ready. The input file is completed with required headers and delete instructions at the end.
Check the following links to see the complete input files:

* BeadsColumnTS.cpp for the Time-Stepping version
* BeadsColumnED.cpp for the Event-Driven
* BeadsPlugin.cpp for the file that contains external plug-in


Results
-------

You can now run in a terminal::

  siconos multiBeadsTS.cpp

and then plot with for example gnuplot::

  gnuplot -persist result.gp

result.gp being a command file (see example in mechanics/MultiBeadsColumn)

Results are given in fig 2, below:

.. figure:: /figures/mechanics/MultiBeads/MultiBeads.*
   :align: center

   fig 2: Result of MultiBeads simulation

.. highlight:: python