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