1#!/usr/bin/env python
2
3from siconos.kernel import SiconosVector
4import sys
5import os
6from siconos.tests_setup import working_dir
7
8v1 = SiconosVector([1.0, 2.0, 3.0])
9
10v2 = SiconosVector()
11
12def test_serialization1():
13    ''' test xml IO'''
14
15    v2.xml_import(v1.xml_export())
16
17    assert(v2.getValue(0) == 1.0)
18    assert(v2.getValue(1) == 2.0)
19    assert(v2.getValue(2) == 3.0)
20
21# question to Maurice: should we keep that or is it completely obsolete?
22#def test_serialization2():
23#    ''' test text IO'''
24#
25#    v2.text_import(v1.text_export())
26#
27#    assert(v2.getValue(0) == 1.0)
28#    assert(v2.getValue(1) == 2.0)
29#    assert(v2.getValue(2) == 3.0)
30#
31
32def test_serialization3():
33    ''' test binary IO'''
34
35    v2.binary_import(v1.binary_export())
36
37    assert(v2.getValue(0) == 1.0)
38    assert(v2.getValue(1) == 2.0)
39    assert(v2.getValue(2) == 3.0)
40
41def test_serialization4():
42    from siconos.kernel import LagrangianLinearTIDS, NewtonImpactNSL, \
43        LagrangianLinearTIR, Interaction, NonSmoothDynamicalSystem, \
44        MoreauJeanOSI, TimeDiscretisation, LCP, TimeStepping
45
46    from numpy import array, eye, empty
47
48    t0 = 0       # start time
49    T = 10       # end time
50    h = 0.005    # time step
51    r = 0.1      # ball radius
52    g = 9.81     # gravity
53    m = 1        # ball mass
54    e = 0.9      # restitution coeficient
55    theta = 0.5  # theta scheme
56
57    #
58    # dynamical system
59    #
60    x = array([1, 0, 0])  # initial position
61    v = array([0, 0, 0])  # initial velocity
62    mass = eye(3)         # mass matrix
63    mass[2, 2] = 3./5 * r * r
64
65    # the dynamical system
66    ball = LagrangianLinearTIDS(x, v, mass)
67
68    # set external forces
69    weight = array([-m * g, 0, 0])
70    ball.setFExtPtr(weight)
71
72    #
73    # Interactions
74    #
75
76    # ball-floor
77    H = array([[1, 0, 0]])
78
79    nslaw = NewtonImpactNSL(e)
80    relation = LagrangianLinearTIR(H)
81    inter = Interaction(nslaw, relation)
82
83    #
84    # Model
85    #
86    first_bouncingBall = NonSmoothDynamicalSystem(t0, T)
87
88    # add the dynamical system to the non smooth dynamical system
89    first_bouncingBall.insertDynamicalSystem(ball)
90
91    # link the interaction and the dynamical system
92    first_bouncingBall.link(inter, ball)
93
94    #
95    # Simulation
96    #
97
98    # (1) OneStepIntegrators
99    OSI = MoreauJeanOSI(theta)
100
101    # (2) Time discretisation --
102    t = TimeDiscretisation(t0, h)
103
104    # (3) one step non smooth problem
105    osnspb = LCP()
106
107    # (4) Simulation setup with (1) (2) (3)
108    s = TimeStepping(first_bouncingBall, t)
109    s.insertIntegrator(OSI)
110    s.insertNonSmoothProblem(osnspb)
111
112    # end of model definition
113
114    #
115    # save and load data from xml and .dat
116    #
117    from siconos.io.io_base import save, load
118    save(s, "bouncingBall.xml")
119
120    s = load("bouncingBall.xml")
121    bouncingBall = s.nonSmoothDynamicalSystem()
122    ball = bouncingBall.dynamicalSystem(ball.number())
123    inter = bouncingBall.interaction(inter.number())
124
125    # the number of time steps
126    N = int((T-t0)/h+1)
127
128    # Get the values to be plotted
129    # ->saved in a matrix dataPlot
130
131    dataPlot = empty((N, 5))
132
133    #
134    # numpy pointers on dense Siconos vectors
135    #
136    q = ball.q()
137    v = ball.velocity()
138    p = ball.p(1)
139    lambda_ = inter.lambda_(1)
140
141    #
142    # initial data
143    #
144    dataPlot[0, 0] = t0
145    dataPlot[0, 1] = q[0]
146    dataPlot[0, 2] = v[0]
147    dataPlot[0, 3] = p[0]
148    dataPlot[0, 4] = lambda_[0]
149
150    k = 1
151
152    # time loop
153    while(s.hasNextEvent()):
154        s.computeOneStep()
155
156        dataPlot[k, 0] = s.nextTime()
157        dataPlot[k, 1] = q[0]
158        dataPlot[k, 2] = v[0]
159        dataPlot[k, 3] = p[0]
160        dataPlot[k, 4] = lambda_[0]
161
162        k += 1
163        print(s.nextTime())
164        s.nextStep()
165
166    #
167    # comparison with the reference file
168    #
169    from siconos.kernel import SimpleMatrix, getMatrix
170    from numpy.linalg import norm
171
172    ref = getMatrix(SimpleMatrix(os.path.join(working_dir,
173                                              "data/result.ref")))
174
175    assert (norm(dataPlot - ref) < 1e-12)
176