GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: include/pinocchio/multibody/data.hxx Lines: 221 227 97.4 %
Date: 2024-01-23 21:41:47 Branches: 326 640 50.9 %

Line Branch Exec Source
1
//
2
// Copyright (c) 2015-2022 CNRS INRIA
3
// Copyright (c) 2015 Wandercraft, 86 rue de Paris 91400 Orsay, France.
4
//
5
6
#ifndef __pinocchio_multibody_data_hxx__
7
#define __pinocchio_multibody_data_hxx__
8
9
#include "pinocchio/spatial/fwd.hpp"
10
#include "pinocchio/multibody/model.hpp"
11
#include "pinocchio/utils/string-generator.hpp"
12
#include "pinocchio/multibody/liegroup/liegroup-algo.hpp"
13
14
/// @cond DEV
15
16
namespace pinocchio
17
{
18
19
  template<typename Scalar, int Options, template<typename,int> class JointCollectionTpl>
20
392
  inline DataTpl<Scalar,Options,JointCollectionTpl>::
21
  DataTpl(const Model & model)
22
  : joints(0)
23
392
  , a((std::size_t)model.njoints,Motion::Zero())
24
392
  , oa((std::size_t)model.njoints,Motion::Zero())
25
392
  , a_gf((std::size_t)model.njoints,Motion::Zero())
26
392
  , oa_gf((std::size_t)model.njoints,Motion::Zero())
27
392
  , v((std::size_t)model.njoints,Motion::Zero())
28
392
  , ov((std::size_t)model.njoints,Motion::Zero())
29
392
  , f((std::size_t)model.njoints,Force::Zero())
30
392
  , of((std::size_t)model.njoints,Force::Zero())
31
392
  , h((std::size_t)model.njoints,Force::Zero())
32
392
  , oh((std::size_t)model.njoints,Force::Zero())
33
392
  , oMi((std::size_t)model.njoints,SE3::Identity())
34
392
  , liMi((std::size_t)model.njoints,SE3::Identity())
35
392
  , tau(VectorXs::Zero(model.nv))
36
392
  , nle(VectorXs::Zero(model.nv))
37
392
  , g(VectorXs::Zero(model.nv))
38
392
  , oMf((std::size_t)model.nframes,SE3::Identity())
39
392
  , Ycrb((std::size_t)model.njoints,Inertia::Zero())
40
392
  , dYcrb((std::size_t)model.njoints,Inertia::Zero())
41
392
  , M(MatrixXs::Zero(model.nv,model.nv))
42
392
  , Minv(MatrixXs::Zero(model.nv,model.nv))
43
392
  , C(MatrixXs::Zero(model.nv,model.nv))
44
392
  , dHdq(Matrix6x::Zero(6,model.nv))
45
392
  , dFdq(Matrix6x::Zero(6,model.nv))
46
392
  , dFdv(Matrix6x::Zero(6,model.nv))
47
392
  , dFda(Matrix6x::Zero(6,model.nv))
48
392
  , SDinv(Matrix6x::Zero(6,model.nv))
49
392
  , UDinv(Matrix6x::Zero(6,model.nv))
50
392
  , IS(MatrixXs::Zero(6,model.nv))
51
392
  , vxI((std::size_t)model.njoints,Inertia::Matrix6::Zero())
52
392
  , Ivx((std::size_t)model.njoints,Inertia::Matrix6::Zero())
53
392
  , B((std::size_t)model.njoints,Inertia::Matrix6::Zero())
54
392
  , oinertias((std::size_t)model.njoints,Inertia::Zero())
55
392
  , oYcrb((std::size_t)model.njoints,Inertia::Zero())
56
392
  , doYcrb((std::size_t)model.njoints,Inertia::Matrix6::Zero())
57
392
  , ddq(VectorXs::Zero(model.nv))
58
392
  , Yaba((std::size_t)model.njoints,Inertia::Matrix6::Zero())
59
392
  , u(VectorXs::Zero(model.nv))
60
392
  , Ag(Matrix6x::Zero(6,model.nv))
61
392
  , dAg(Matrix6x::Zero(6,model.nv))
62
  , hg(Force::Zero())
63
  , dhg(Force::Zero())
64
  , Ig(Inertia::Zero())
65
784
  , Fcrb((std::size_t)model.njoints,Matrix6x::Zero(6,model.nv))
66
  , lastChild((std::size_t)model.njoints,-1)
67
  , nvSubtree((std::size_t)model.njoints,-1)
68
  , start_idx_v_fromRow((std::size_t)model.nv,-1)
69
  , end_idx_v_fromRow((std::size_t)model.nv,-1)
70
392
  , U(MatrixXs::Identity(model.nv,model.nv))
71
392
  , D(VectorXs::Zero(model.nv))
72
392
  , Dinv(VectorXs::Zero(model.nv))
73
392
  , tmp(VectorXs::Zero(model.nv))
74
  , parents_fromRow((std::size_t)model.nv,-1)
75
392
  , supports_fromRow((std::size_t)model.nv)
76
  , nvSubtree_fromRow((std::size_t)model.nv,-1)
77
392
  , J(Matrix6x::Zero(6,model.nv))
78
392
  , dJ(Matrix6x::Zero(6,model.nv))
79
392
  , ddJ(Matrix6x::Zero(6, model.nv))
80
392
  , psid(Matrix6x::Zero(6, model.nv))
81
392
  , psidd(Matrix6x::Zero(6, model.nv))
82
392
  , dVdq(Matrix6x::Zero(6,model.nv))
83
392
  , dAdq(Matrix6x::Zero(6,model.nv))
84
392
  , dAdv(Matrix6x::Zero(6,model.nv))
85
392
  , dtau_dq(MatrixXs::Zero(model.nv,model.nv))
86
392
  , dtau_dv(MatrixXs::Zero(model.nv,model.nv))
87
392
  , ddq_dq(MatrixXs::Zero(model.nv,model.nv))
88
392
  , ddq_dv(MatrixXs::Zero(model.nv,model.nv))
89
392
  , iMf((std::size_t)model.njoints,SE3::Identity())
90
392
  , com((std::size_t)model.njoints,Vector3::Zero())
91
392
  , vcom((std::size_t)model.njoints,Vector3::Zero())
92
392
  , acom((std::size_t)model.njoints,Vector3::Zero())
93
392
  , mass((std::size_t)model.njoints,(Scalar)(-1))
94
392
  , Jcom(Matrix3x::Zero(3,model.nv))
95
  , kinetic_energy((Scalar)-1)
96
  , potential_energy((Scalar)-1)
97
  , JMinvJt()
98
  , llt_JMinvJt()
99
  , lambda_c()
100
392
  , sDUiJt(MatrixXs::Zero(model.nv,model.nv))
101
392
  , torque_residual(VectorXs::Zero(model.nv))
102
392
  , dq_after(VectorXs::Zero(model.nv))
103
  , impulse_c()
104
392
  , staticRegressor(Matrix3x::Zero(3,4*(model.njoints-1)))
105
  , bodyRegressor(BodyRegressorType::Zero())
106
392
  , jointTorqueRegressor(MatrixXs::Zero(model.nv,10*(model.njoints-1)))
107
#if EIGEN_VERSION_AT_LEAST(3,2,90) && !EIGEN_VERSION_AT_LEAST(3,2,93)
108
  , kinematic_hessians(6,std::max(1,model.nv),std::max(1,model.nv)) // the minimum size should be 1 for compatibility reasons
109
  , d2tau_dqdq(std::max(1,model.nv),std::max(1,model.nv),std::max(1,model.nv)) // the minimum size should be 1 for compatibility reasons
110
  , d2tau_dvdv(std::max(1,model.nv),std::max(1,model.nv),std::max(1,model.nv)) // the minimum size should be 1 for compatibility reasons
111
  , d2tau_dqdv(std::max(1,model.nv),std::max(1,model.nv),std::max(1,model.nv)) // the minimum size should be 1 for compatibility reasons
112
  , d2tau_dadq(std::max(1,model.nv),std::max(1,model.nv),std::max(1,model.nv)) // the minimum size should be 1 for compatibility reasons
113
 #else
114
392
   , kinematic_hessians(6,model.nv,model.nv)
115
392
   , d2tau_dqdq(model.nv,model.nv,model.nv)
116
392
   , d2tau_dvdv(model.nv,model.nv,model.nv)
117
392
   , d2tau_dqdv(model.nv,model.nv,model.nv)
118




















































































392
   , d2tau_dadq(model.nv,model.nv,model.nv)
119
#endif
120
  {
121
    typedef typename Model::JointIndex JointIndex;
122
123
    /* Create data structure associated to the joints */
124
10271
    for(JointIndex i=0;i<(JointIndex)(model.njoints);++i)
125

9879
      joints.push_back(CreateJointData<Scalar,Options,JointCollectionTpl>::run(model.joints[i]));
126
127
    /* Init for CRBA */
128

392
    M.setZero(); Minv.setZero();
129
10271
    for(JointIndex i=0;i<(JointIndex)(model.njoints);++i)
130
9879
    { Fcrb[i].resize(6,model.nv); }
131
132
392
    computeLastChild(model);
133
134
    /* Init for Coriolis */
135
136
    /* Init for Cholesky */
137
392
    computeParents_fromRow(model);
138
392
    computeSupports_fromRow(model);
139
140
    /* Init universe states relatively to itself */
141
392
    a_gf[0] = -model.gravity;
142
143
392
    kinematic_hessians.setZero();
144
392
    d2tau_dqdq.setZero();
145
392
    d2tau_dvdv.setZero();
146
392
    d2tau_dqdv.setZero();
147
392
    d2tau_dadq.setZero();
148
149
392
  }
150
151
  template<typename Scalar, int Options, template<typename,int> class JointCollectionTpl>
152
392
  inline void DataTpl<Scalar,Options,JointCollectionTpl>::
153
  computeLastChild(const Model & model)
154
  {
155
    typedef typename Model::Index Index;
156
157
392
    std::fill(lastChild.begin(),lastChild.end(),-1);
158
10271
    for( int i=model.njoints-1;i>=0;--i )
159
    {
160
9879
      if(lastChild[(Index)i] == -1) lastChild[(Index)i] = i;
161
9879
      const Index & parent = model.parents[(Index)i];
162
9879
      lastChild[parent] = std::max<int>(lastChild[(Index)i],lastChild[parent]);
163
164
19758
      nvSubtree[(Index)i]
165
9879
      = idx_v(model.joints[(Index)lastChild[(Index)i]]) + nv(model.joints[(Index)lastChild[(Index)i]])
166
9879
      - idx_v(model.joints[(Index)i]);
167
    }
168
392
  }
169
170
  template<typename Scalar, int Options, template<typename,int> class JointCollectionTpl>
171
392
  inline void DataTpl<Scalar,Options,JointCollectionTpl>
172
  ::computeParents_fromRow(const Model & model)
173
  {
174
    typedef typename Model::Index Index;
175
176
9879
    for(Index joint=1;joint<(Index)(model.njoints);joint++)
177
    {
178
9487
      const Index & parent = model.parents[joint];
179
9487
      const int nvj    = nv   (model.joints[joint]);
180
9487
      const int idx_vj = idx_v(model.joints[joint]);
181
182

9487
      assert(idx_vj >= 0 && idx_vj < model.nv);
183
184
9487
      if(parent>0) parents_fromRow[(size_t)idx_vj] = idx_v(model.joints[parent])+nv(model.joints[parent])-1;
185
392
      else         parents_fromRow[(size_t)idx_vj] = -1;
186
9487
      nvSubtree_fromRow[(size_t)idx_vj] = nvSubtree[joint];
187
188
9487
      start_idx_v_fromRow[(size_t)idx_vj] = idx_vj;
189
9487
      end_idx_v_fromRow[(size_t)idx_vj] = idx_vj+nvj-1;
190
11299
      for(int row=1;row<nvj;++row)
191
      {
192
1812
        parents_fromRow[(size_t)(idx_vj+row)] = idx_vj+row-1;
193
1812
        nvSubtree_fromRow[(size_t)(idx_vj+row)] = nvSubtree[joint]-row;
194
1812
        start_idx_v_fromRow[(size_t)(idx_vj+row)] = start_idx_v_fromRow[(size_t)idx_vj];
195
1812
        end_idx_v_fromRow[(size_t)(idx_vj+row)] = end_idx_v_fromRow[(size_t)idx_vj];
196
      }
197
    }
198
392
  }
199
200
  template<typename Scalar, int Options, template<typename,int> class JointCollectionTpl>
201
392
  inline void DataTpl<Scalar,Options,JointCollectionTpl>
202
  ::computeSupports_fromRow(const Model & model)
203
  {
204
    typedef typename Model::JointIndex JointIndex;
205
206
9879
    for(JointIndex joint_id = 1;
207
9879
        joint_id < (JointIndex)(model.njoints);
208
        joint_id++)
209
    {
210
9487
      const int nvj    = nv   (model.joints[joint_id]);
211
9487
      const int idx_vj = idx_v(model.joints[joint_id]);
212
213

9487
      assert(idx_vj >= 0 && idx_vj < model.nv);
214
215
9487
      const int parent_fromRow = parents_fromRow[(size_t)idx_vj];
216
217
9487
      if(parent_fromRow >= 0)
218
9093
        supports_fromRow[(size_t)idx_vj] = supports_fromRow[(size_t)parent_fromRow];
219
220
9487
      supports_fromRow[(size_t)idx_vj].push_back(idx_vj);
221
222
11299
      for(int row = 1; row < nvj; ++row)
223
      {
224
1812
        supports_fromRow[(size_t)(idx_vj+row)] = supports_fromRow[(size_t)(idx_vj+row-1)];
225
1812
        supports_fromRow[(size_t)(idx_vj+row)].push_back(idx_vj+row);
226
      }
227
    }
228
392
  }
229
230
  template<typename Scalar, int Options, template<typename,int> class JointCollectionTpl>
231
10
  bool operator==(const DataTpl<Scalar,Options,JointCollectionTpl> & data1,
232
                  const DataTpl<Scalar,Options,JointCollectionTpl> & data2)
233
  {
234
10
    bool value =
235
10
       data1.joints == data2.joints
236
10
    && data1.a == data2.a
237
10
    && data1.oa == data2.oa
238
10
    && data1.a_gf == data2.a_gf
239
10
    && data1.oa_gf == data2.oa_gf
240
10
    && data1.v == data2.v
241
10
    && data1.ov == data2.ov
242
10
    && data1.f == data2.f
243
10
    && data1.of == data2.of
244
10
    && data1.h == data2.h
245
10
    && data1.oh == data2.oh
246
10
    && data1.oMi == data2.oMi
247
9
    && data1.liMi == data2.liMi
248
9
    && data1.tau == data2.tau
249
9
    && data1.nle == data2.nle
250
9
    && data1.g == data2.g
251
9
    && data1.oMf == data2.oMf
252
9
    && data1.Ycrb == data2.Ycrb
253
9
    && data1.dYcrb == data2.dYcrb
254
9
    && data1.M == data2.M
255
9
    && data1.Minv == data2.Minv
256
9
    && data1.C == data2.C
257
9
    && data1.dHdq == data2.dHdq
258
9
    && data1.dFdq == data2.dFdq
259
9
    && data1.dFdv == data2.dFdv
260
9
    && data1.dFda == data2.dFda
261
9
    && data1.SDinv == data2.SDinv
262
9
    && data1.UDinv == data2.UDinv
263
9
    && data1.IS == data2.IS
264
9
    && data1.vxI == data2.vxI
265
9
    && data1.Ivx == data2.Ivx
266
9
    && data1.oinertias == data2.oinertias
267
9
    && data1.oYcrb == data2.oYcrb
268
9
    && data1.doYcrb == data2.doYcrb
269
9
    && data1.ddq == data2.ddq
270
9
    && data1.Yaba == data2.Yaba
271
9
    && data1.u == data2.u
272
9
    && data1.Ag == data2.Ag
273
9
    && data1.dAg == data2.dAg
274
9
    && data1.hg == data2.hg
275
9
    && data1.dhg == data2.dhg
276
9
    && data1.Ig == data2.Ig
277
9
    && data1.Fcrb == data2.Fcrb
278
9
    && data1.lastChild == data2.lastChild
279
9
    && data1.nvSubtree == data2.nvSubtree
280
9
    && data1.start_idx_v_fromRow == data2.start_idx_v_fromRow
281
9
    && data1.end_idx_v_fromRow == data2.end_idx_v_fromRow
282
9
    && data1.U == data2.U
283
9
    && data1.D == data2.D
284
9
    && data1.Dinv == data2.Dinv
285
9
    && data1.parents_fromRow == data2.parents_fromRow
286
9
    && data1.supports_fromRow == data2.supports_fromRow
287
9
    && data1.nvSubtree_fromRow == data2.nvSubtree_fromRow
288
9
    && data1.J == data2.J
289
9
    && data1.dJ == data2.dJ
290
9
    && data1.ddJ == data2.ddJ
291
9
    && data1.psid == data2.psid
292
9
    && data1.psidd == data2.psidd
293
9
    && data1.dVdq == data2.dVdq
294
9
    && data1.dAdq == data2.dAdq
295
9
    && data1.dAdv == data2.dAdv
296
9
    && data1.dtau_dq == data2.dtau_dq
297
9
    && data1.dtau_dv == data2.dtau_dv
298
9
    && data1.ddq_dq == data2.ddq_dq
299
9
    && data1.ddq_dv == data2.ddq_dv
300
9
    && data1.iMf == data2.iMf
301
9
    && data1.com == data2.com
302
9
    && data1.vcom == data2.vcom
303
9
    && data1.acom == data2.acom
304
9
    && data1.mass == data2.mass
305
9
    && data1.Jcom == data2.Jcom
306
9
    && data1.kinetic_energy == data2.kinetic_energy
307
9
    && data1.potential_energy == data2.potential_energy
308
9
    && data1.JMinvJt == data2.JMinvJt
309
9
    && data1.lambda_c == data2.lambda_c
310
9
    && data1.torque_residual == data2.torque_residual
311
9
    && data1.dq_after == data2.dq_after
312
9
    && data1.impulse_c == data2.impulse_c
313
9
    && data1.staticRegressor == data2.staticRegressor
314
9
    && data1.bodyRegressor == data2.bodyRegressor
315

20
    && data1.jointTorqueRegressor == data2.jointTorqueRegressor
316
    ;
317
318
    // operator== for Eigen::Tensor provides an Expression which might be not evaluated as a boolean
319

20
    value &= Tensor<bool, 0>((data1.kinematic_hessians == data2.kinematic_hessians).all())(0)
320



20
    && Tensor<bool, 0>((data1.d2tau_dqdq == data2.d2tau_dqdq).all())(0)
321



20
    && Tensor<bool, 0>((data1.d2tau_dvdv == data2.d2tau_dvdv).all())(0)
322



20
    && Tensor<bool, 0>((data1.d2tau_dqdv == data2.d2tau_dqdv).all())(0)
323





30
    && Tensor<bool, 0>((data1.d2tau_dadq == data2.d2tau_dadq).all())(0);
324
325
10
    return value;
326
  }
327
328
  template<typename Scalar, int Options, template<typename,int> class JointCollectionTpl>
329
1
  bool operator!=(const DataTpl<Scalar,Options,JointCollectionTpl> & data1,
330
                  const DataTpl<Scalar,Options,JointCollectionTpl> & data2)
331
  {
332
1
    return !(data1 == data2);
333
  }
334
335
} // namespace pinocchio
336
337
/// @endcond
338
339
#endif // ifndef __pinocchio_multibody_data_hxx__