1 #include "SiconosConfig.h"
2 #define WITH_SERIALIZATION
3
4 #ifdef HAVE_SICONOS_MECHANICS
5 #include "MechanicsIO.hpp"
6 #endif
7
8 #include "KernelTest.hpp"
9 #include "SiconosKernel.hpp"
10
11 #include <boost/numeric/bindings/ublas/matrix.hpp>
12 #include <boost/numeric/bindings/ublas/vector.hpp>
13 #include <boost/numeric/bindings/ublas/vector_sparse.hpp>
14 #include <boost/numeric/bindings/ublas/matrix_sparse.hpp>
15
16 #define DEBUG_MESSAGES 1
17 #include "SiconosFull.hpp"
18
19 #include <boost/archive/binary_iarchive.hpp>
20 #include <boost/archive/binary_oarchive.hpp>
21
22 #include <boost/archive/xml_iarchive.hpp>
23 #include <boost/archive/xml_oarchive.hpp>
24
25 CPPUNIT_TEST_SUITE_REGISTRATION(KernelTest);
26
setUp()27 void KernelTest::setUp()
28 {
29 BBxml = "BouncingBall1.xml";
30 }
31
tearDown()32 void KernelTest::tearDown() {};
33
t0()34 void KernelTest::t0()
35 {
36 SP::SiconosVector q(new SiconosVector(3));
37 SP::SiconosVector q0(new SiconosVector(3));
38
39 (*q)(0) = 1.0;
40 (*q)(1) = 2.0;
41 (*q)(2) = 2.0;
42
43
44 std::ofstream ofs("Kernelt0.xml");
45 {
46 boost::archive::xml_oarchive oa(ofs);
47 oa.register_type(static_cast<SiconosVector*>(nullptr));
48 oa << NVP(q);
49 }
50
51 std::ifstream ifs("Kernelt0.xml");
52 {
53 boost::archive::xml_iarchive ia(ifs);
54 ia.register_type(static_cast<SiconosVector*>(nullptr));
55 ia >> NVP(q0);
56 }
57
58 CPPUNIT_ASSERT(*q0 == *q);
59 }
60
61
t1()62 void KernelTest::t1()
63 {
64 SP::SiconosMatrix m1(new SimpleMatrix(3, 3));
65 SP::SimpleMatrix m2(new SimpleMatrix(3, 3));
66
67 m1->eye();
68 (*m1)(1, 0) = 3.0;
69 (*m1)(2, 1) = -7;
70
71 std::ofstream ofs("Kernelt1.xml");
72 {
73 boost::archive::xml_oarchive oa(ofs);
74 oa.register_type(static_cast<SimpleMatrix*>(nullptr));
75 oa << NVP(m1);
76 }
77
78 std::ifstream ifs("Kernelt1.xml");
79 {
80 boost::archive::xml_iarchive ia(ifs);
81 ia.register_type(static_cast<SimpleMatrix*>(nullptr));
82 ia >> NVP(m2);
83 }
84
85 m1->display();
86 m2->display();
87
88 CPPUNIT_ASSERT(*m1 == *m2);
89
90 }
91
t2()92 void KernelTest::t2()
93 {
94 SP::SiconosMatrix m(new SimpleMatrix(3, 3));
95 SP::SiconosVector v(new SiconosVector(3));
96 SP::SiconosVector q(new SiconosVector(3));
97
98 m->eye();
99
100
101 SP::DynamicalSystem ds1(new LagrangianDS(q, v, m));
102 SP::DynamicalSystem ds2(new LagrangianDS(q, v, m));
103
104 std::ofstream ofs("Kernelt2.xml");
105 {
106 boost::archive::xml_oarchive oa(ofs);
107 oa.register_type(static_cast<SimpleMatrix*>(nullptr));
108 oa.register_type(static_cast<SiconosVector*>(nullptr));
109 oa.register_type(static_cast<LagrangianDS*>(nullptr));
110 oa << NVP(ds1);
111 }
112
113 std::ifstream ifs("Kernelt2.xml");
114 {
115 boost::archive::xml_iarchive ia(ifs);
116 ia.register_type(static_cast<SimpleMatrix*>(nullptr));
117 ia.register_type(static_cast<SiconosVector*>(nullptr));
118 ia.register_type(static_cast<LagrangianDS*>(nullptr));
119 ia >> NVP(ds2);
120 }
121
122 CPPUNIT_ASSERT(*(std::static_pointer_cast<LagrangianDS>(ds1)->mass())
123 == *(std::static_pointer_cast<LagrangianDS>(ds2)->mass()));
124 CPPUNIT_ASSERT(*(std::static_pointer_cast<LagrangianDS>(ds1)->q())
125 == *(std::static_pointer_cast<LagrangianDS>(ds2)->q()));
126 CPPUNIT_ASSERT(*(std::static_pointer_cast<LagrangianDS>(ds1)->velocity())
127 == *(std::static_pointer_cast<LagrangianDS>(ds2)->velocity()));
128
129 }
130
131
t3()132 void KernelTest::t3()
133 {
134 SP::SolverOptions so(solver_options_create(SICONOS_FRICTION_3D_NSN_AC),
135 solver_options_delete);
136 SP::SolverOptions sor(new SolverOptions);
137 so->numberOfInternalSolvers = 1;
138 so->internalSolvers = (SolverOptions**) calloc(1, sizeof(SolverOptions*));
139 so->internalSolvers[0] = solver_options_create(SICONOS_FRICTION_3D_NSN_AC);
140
141 std::ofstream ofs("SolverOptions.xml");
142 {
143 boost::archive::xml_oarchive oa(ofs);
144 oa << NVP(so);
145 }
146
147 std::ifstream ifs("SolverOptions.xml");
148 {
149 boost::archive::xml_iarchive ia(ifs);
150 ia >> NVP(sor);
151 }
152
153 CPPUNIT_ASSERT((so->iSize == sor->iSize));
154 }
155
156 // void KernelTest::t4()
157 // {
158 // SP::SiconosMatrix m(new SimpleMatrix(3, 3));
159 // SP::SiconosVector v(new SiconosVector(3));
160 // SP::SiconosVector q(new SiconosVector(3));
161
162 // m->eye();
163
164
165 // SP::DynamicalSystem ds1(new LagrangianDS(q, v, m));
166 // SP::DynamicalSystem ds2(new LagrangianDS(q, v, m));
167
168 // SP::DynamicalSystemsSet dsset(new DynamicalSystemsSet());
169
170 // dsset->insert(ds1);
171 // dsset->insert(ds2);
172
173 // std::ofstream ofs("t4.xml");
174 // {
175 // boost::archive::xml_oarchive oa(ofs);
176 // oa.register_type(static_cast<SimpleMatrix*>(nullptr));
177 // oa.register_type(static_cast<SiconosVector*>(nullptr));
178 // oa.register_type(static_cast<LagrangianDS*>(nullptr));
179 // oa << NVP(dsset);
180 // }
181
182 // SP::DynamicalSystemsSet dssetfromfile(new DynamicalSystemsSet());
183
184 // std::ifstream ifs("t4.xml");
185 // {
186 // boost::archive::xml_iarchive ia(ifs);
187 // ia.register_type(static_cast<SimpleMatrix*>(nullptr));
188 // ia.register_type(static_cast<SiconosVector*>(nullptr));
189 // ia.register_type(static_cast<LagrangianDS*>(nullptr));
190 // ia >> NVP(dssetfromfile);
191 // }
192
193 // }
194
195
196 #include "SiconosRestart.hpp"
197
198 using namespace std;
t5()199 void KernelTest::t5()
200 {
201
202 // ================= Creation of the model =======================
203
204 // User-defined main parameters
205 unsigned int nDof = 3; // degrees of freedom for the ball
206 double t0 = 0; // initial computation time
207 double T = 10; // final computation time
208 double h = 0.005; // time step
209 double position_init = 1.0; // initial position for lowest bead.
210 double velocity_init = 0.0; // initial velocity for lowest bead.
211 double theta = 0.5; // theta for MoreauJeanOSI integrator
212 double R = 0.1; // Ball radius
213 double m = 1; // Ball mass
214 double g = 9.81; // Gravity
215 // -------------------------
216 // --- Dynamical systems ---
217 // -------------------------
218
219 cout << "====> Model loading ..." << endl << endl;
220
221 SP::SiconosMatrix Mass(new SimpleMatrix(nDof, nDof));
222 (*Mass)(0, 0) = m;
223 (*Mass)(1, 1) = m;
224 (*Mass)(2, 2) = 3. / 5 * m * R * R;
225
226 // -- Initial positions and velocities --
227 SP::SiconosVector q0(new SiconosVector(nDof));
228 SP::SiconosVector v0(new SiconosVector(nDof));
229 (*q0)(0) = position_init;
230 (*v0)(0) = velocity_init;
231
232 // -- The dynamical system --
233 SP::LagrangianLinearTIDS ball(new LagrangianLinearTIDS(q0, v0, Mass));
234
235 // -- Set external forces (weight) --
236 SP::SiconosVector weight(new SiconosVector(nDof));
237 (*weight)(0) = -m * g;
238 ball->setFExtPtr(weight);
239
240 // --------------------
241 // --- Interactions ---
242 // --------------------
243
244 // -- nslaw --
245 double e = 0.9;
246
247 // Interaction ball-floor
248 //
249 SP::SimpleMatrix H(new SimpleMatrix(1, nDof));
250 (*H)(0, 0) = 1.0;
251
252 SP::NonSmoothLaw nslaw(new NewtonImpactNSL(e));
253 SP::Relation relation(new LagrangianLinearTIR(H));
254
255 SP::Interaction inter(new Interaction(nslaw, relation));
256
257 // -------------
258 // --- Model ---
259 // -------------
260 SP::NonSmoothDynamicalSystem bouncingBall(new NonSmoothDynamicalSystem(t0, T));
261
262 // add the dynamical system in the non smooth dynamical system
263 bouncingBall->insertDynamicalSystem(ball);
264
265 // link the interaction and the dynamical system
266 bouncingBall->link(inter, ball);
267
268
269 // ------------------
270 // --- Simulation ---
271 // ------------------
272
273 // -- (1) OneStepIntegrators --
274 SP::MoreauJeanOSI OSI(new MoreauJeanOSI(theta));
275
276 // -- (2) Time discretisation --
277 SP::TimeDiscretisation t(new TimeDiscretisation(t0, h));
278
279 // -- (3) one step non smooth problem
280 SP::OneStepNSProblem osnspb(new LCP());
281
282 // -- (4) Simulation setup with (1) (2) (3)
283 SP::TimeStepping s(new TimeStepping(bouncingBall, t, OSI, osnspb));
284 s->associate(OSI, ball);
285
286 // =========================== End of model definition ===========================
287
288 // ================================= Computation =================================
289
290 Siconos::save(s, BBxml);
291
292 SP::Simulation simFromFile = Siconos::load(BBxml);
293 SP::NonSmoothDynamicalSystem bouncingBallFromFile =
294 simFromFile->nonSmoothDynamicalSystem();
295
296 CPPUNIT_ASSERT((bouncingBallFromFile->t0() == bouncingBall->t0()));
297 // in depth comparison?
298
299 // now we should try to run the bouncing ball from file
300
301 // BUT: non serialized members => must be initialized or serialized
302
303 }
304
305
t6()306 void KernelTest::t6()
307 {
308 SP::Simulation sim = Siconos::load(BBxml);
309
310 try
311 {
312 SP::NonSmoothDynamicalSystem bouncingBall = sim->nonSmoothDynamicalSystem();
313
314 double T = bouncingBall->finalT();
315 double t0 = bouncingBall->t0();
316 double h = sim->timeStep();
317 int N = (int)((T - t0) / h); // Number of time steps
318
319 SP::DynamicalSystemsGraph dsg =
320 bouncingBall->topology()->dSG(0);
321
322 SP::LagrangianDS ball = std::static_pointer_cast<LagrangianDS>
323 (dsg->bundle(*(dsg->begin())));
324
325 SP::TimeStepping s = std::static_pointer_cast<TimeStepping>(sim);
326 SP::Interaction inter;
327 InteractionsGraph::VIterator ui, uiend;
328 SP::InteractionsGraph indexSet0 = bouncingBall->topology()->indexSet(0);
329 for(std::tie(ui, uiend) = indexSet0->vertices(); ui != uiend; ++ui)
330 inter = indexSet0->bundle(*ui);
331
332
333 // --- Get the values to be plotted ---
334 // -> saved in a matrix dataPlot
335 unsigned int outputSize = 5;
336 SimpleMatrix dataPlot(N + 1, outputSize);
337
338
339
340 SP::SiconosVector q = ball->q();
341 SP::SiconosVector v = ball->velocity();
342 SP::SiconosVector p = ball->p(1);
343 SP::SiconosVector lambda = inter->lambda(1);
344
345 dataPlot(0, 0) = bouncingBall->t0();
346 dataPlot(0, 1) = (*q)(0);
347 dataPlot(0, 2) = (*v)(0);
348 dataPlot(0, 3) = (*p)(0);
349 dataPlot(0, 4) = (*lambda)(0);
350 // --- Time loop ---
351 cout << "====> Start computation ... " << endl << endl;
352 // ==== Simulation loop - Writing without explicit event handling =====
353 int k = 1;
354 while(s->hasNextEvent())
355 {
356 s->computeOneStep();
357
358 // --- Get values to be plotted ---
359 dataPlot(k, 0) = s->nextTime();
360 dataPlot(k, 1) = (*q)(0);
361 dataPlot(k, 2) = (*v)(0);
362 dataPlot(k, 3) = (*p)(0);
363 dataPlot(k, 4) = (*lambda)(0);
364 s->nextStep();
365 k++;
366 }
367 cout << endl << "End of computation - Number of iterations done: " << k - 1 << endl;
368
369 // --- Output files ---
370 cout << "====> Output file writing ..." << endl;
371 dataPlot.resize(k, outputSize);
372 ioMatrix::write("result.dat", "ascii", dataPlot, "noDim");
373 // Comparison with a reference file
374 SimpleMatrix dataPlotRef(dataPlot);
375 dataPlotRef.zero();
376 ioMatrix::read("result.ref", "ascii", dataPlotRef);
377
378 if((dataPlot - dataPlotRef).normInf() > 1e-12)
379 {
380 std::cout <<
381 "Warning. The results is rather different from the reference file :"
382 <<
383 (dataPlot - dataPlotRef).normInf()
384 <<
385 std::endl;
386 CPPUNIT_ASSERT(false);
387 }
388
389 }
390
391 catch(...)
392 {
393 Siconos::exception::process();
394 CPPUNIT_ASSERT(false);
395 }
396 }
397
398 #ifdef HAVE_SICONOS_MECHANICS
399 #include <Disk.hpp>
400
t7()401 void KernelTest::t7()
402 {
403
404 SP::DynamicalSystem ds1, ds2;
405
406 // Must be size=1, cannot deserialize a LagrangianDS with _ndof==0
407 SP::SiconosVector q(new SiconosVector(1));
408 SP::SiconosVector v(new SiconosVector(1));
409
410 ds1.reset(new Disk(1, 1, q, v));
411
412 ds2.reset(new Disk(2, 2, q, v));
413
414 std::ofstream ofs("Kernelt7.xml");
415 {
416 boost::archive::xml_oarchive oa(ofs);
417 siconos_io_register_Mechanics(oa);
418 oa << NVP(ds1);
419 }
420
421 std::ifstream ifs("Kernelt7.xml");
422 {
423 boost::archive::xml_iarchive ia(ifs);
424 siconos_io_register_Mechanics(ia);
425 ia >> NVP(ds2);
426 }
427
428 CPPUNIT_ASSERT(std::static_pointer_cast<Disk>(ds1)->getRadius() ==
429 std::static_pointer_cast<Disk>(ds2)->getRadius());
430 }
431
t8()432 void KernelTest::t8()
433 {
434 SP::DynamicalSystem ds1, ds2;
435
436 SP::SiconosVector q(new SiconosVector(3));
437 SP::SiconosVector v(new SiconosVector(3));
438
439 (*q)(0) = 0.;
440 (*q)(1) = 1.;
441 (*q)(2) = 1.;
442
443 (*v)(0) = 0;
444 (*v)(1) = 0;
445 (*v)(2) = 10.;
446
447 SP::NonSmoothDynamicalSystem nsds(new NonSmoothDynamicalSystem(0,10));
448
449 ds1.reset(new Disk(1, 1, q, v));
450 ds2.reset(new Disk(2, 2, q, v));
451
452 nsds->insertDynamicalSystem(ds1);
453 nsds->insertDynamicalSystem(ds2);
454
455 MechanicsIO IO;
456
457 SP::SimpleMatrix positions = IO.positions(*nsds);
458 SP::SimpleMatrix velocities = IO.velocities(*nsds);
459
460 //ids
461 CPPUNIT_ASSERT((*positions)(0,0) == 0);
462 CPPUNIT_ASSERT((*velocities)(0,0) == 0);
463 CPPUNIT_ASSERT((*positions)(0,0) == ds1->number());
464 CPPUNIT_ASSERT((*velocities)(0,0) == ds1->number());
465
466 CPPUNIT_ASSERT((*positions)(1,0) == 1);
467 CPPUNIT_ASSERT((*velocities)(1,0) == 1);
468 CPPUNIT_ASSERT((*positions)(1,0) == ds2->number());
469 CPPUNIT_ASSERT((*velocities)(1,0) == ds2->number());
470
471 CPPUNIT_ASSERT((*positions)(0,1) == 0.);
472 CPPUNIT_ASSERT((*velocities)(0,1) == 0.);
473 CPPUNIT_ASSERT((*positions)(0,2) == 1.);
474 CPPUNIT_ASSERT((*positions)(1,2) == 1.);
475 CPPUNIT_ASSERT((*velocities)(0,3) == 10.);
476 CPPUNIT_ASSERT((*velocities)(1,3) == 10.);
477
478 }
479 #endif
480
t9()481 void KernelTest::t9()
482 {
483 try
484 {
485 // Serialize and deserialize an NSDS with T=inf
486 // (a possible failure case for xml archives)
487 double t0 = 0.0;
488 double T = std::numeric_limits<double>::infinity();
489 SP::NonSmoothDynamicalSystem nsds1(new NonSmoothDynamicalSystem(t0, T));
490 SP::NonSmoothDynamicalSystem nsds2;
491
492 std::ofstream ofs("Kernelt9.xml");
493 {
494 boost::archive::xml_oarchive oa(ofs);
495 siconos_io_register_Kernel(oa);
496 oa << NVP(nsds1);
497 }
498
499 std::ifstream ifs("Kernelt9.xml");
500 {
501 boost::archive::xml_iarchive ia(ifs);
502 siconos_io_register_Kernel(ia);
503 ia >> NVP(nsds2);
504 }
505
506 CPPUNIT_ASSERT(nsds1->finalT() == nsds2->finalT());
507 }
508 catch(...)
509 {
510 Siconos::exception::process();
511 CPPUNIT_ASSERT(false);
512 }
513 }
514