1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Alessandro Tasora, Radu Serban
13 // =============================================================================
14
15 #ifndef CHLOAD_H
16 #define CHLOAD_H
17
18 #include "chrono/physics/ChLoader.h"
19 #include "chrono/physics/ChLoaderU.h"
20 #include "chrono/physics/ChLoaderUV.h"
21 #include "chrono/physics/ChLoaderUVW.h"
22 #include "chrono/physics/ChObject.h"
23 #include "chrono/solver/ChKblockGeneric.h"
24 #include "chrono/solver/ChSystemDescriptor.h"
25 #include "chrono/timestepper/ChState.h"
26
27 namespace chrono {
28
29 /// Utility class for storing jacobian matrices.
30 /// This is automatically managed by the ChLoad, if needed
31 /// (ie. for stiff loads)
32
33 class ChApi ChLoadJacobians {
34 public:
35 ChKblockGeneric KRM; ///< sum of K,R,M, with pointers to sparse variables
36 ChMatrixDynamic<double> K; ///< dQ/dx
37 ChMatrixDynamic<double> R; ///< dQ/dv
38 ChMatrixDynamic<double> M; ///< dQ/da
39
40 /// Set references to the constrained objects, each of ChVariables type,
41 /// automatically creating/resizing K matrix if needed.
42 void SetVariables(std::vector<ChVariables*> mvariables);
43 };
44
45 // -----------------------------------------------------------------------------
46
47 /// Base class for loads.
48 /// This class can be inherited to implement applied loads to bodies,
49 /// finite elements, etc.
50 /// A load is an object that might wrap one or more ChLoader objects, whose
51 /// value is dependent.
52 /// It implements functionalities to perform automatic differentiation of
53 /// the load so it optionally can compute the jacobian (the tangent stiffness
54 /// matrix of the load) that can be used in implicit integrators, statics, etc.
55
56 class ChApi ChLoadBase : public ChObj {
57 protected:
58 ChLoadJacobians* jacobians;
59
60 public:
61 ChLoadBase();
62
63 virtual ~ChLoadBase();
64
65 /// Gets the number of DOFs affected by this load (position part)
66 virtual int LoadGet_ndof_x() = 0;
67
68 /// Gets the number of DOFs affected by this load (speed part)
69 virtual int LoadGet_ndof_w() = 0;
70
71 /// Gets all the current DOFs packed in a single vector (position part)
72 virtual void LoadGetStateBlock_x(ChState& mD) = 0;
73
74 /// Gets all the current DOFs packed in a single vector (speed part)
75 virtual void LoadGetStateBlock_w(ChStateDelta& mD) = 0;
76
77 /// Increment a packed state (ex. as obtained by LoadGetStateBlock_x()) by a given packed state-delta.
78 /// Compute: x_new = x + dw. Ex. this is called by the BDF numerical differentiation routine that computes jacobian
79 /// in the default ComputeJacobian() fallback, if not overriding ComputeJacobian() with an analytical form.
80 virtual void LoadStateIncrement(const ChState& x, const ChStateDelta& dw, ChState& x_new) = 0;
81
82 /// Number of coordinates in the interpolated field, ex=3 for a
83 /// tetrahedron finite element or a cable, = 1 for a thermal problem, etc.
84 virtual int LoadGet_field_ncoords() = 0;
85
86 /// Compute Q, the generalized load(s).
87 /// Where Q is stored depends on children classes.
88 /// Called automatically at each Update().
89 virtual void ComputeQ(ChState* state_x, ///< state position to evaluate Q
90 ChStateDelta* state_w ///< state speed to evaluate Q
91 ) = 0;
92
93 /// Compute the K=-dQ/dx, R=-dQ/dv , M=-dQ/da jacobians.
94 /// Called automatically at each Update().
95 virtual void ComputeJacobian(ChState* state_x, ///< state position to evaluate jacobians
96 ChStateDelta* state_w, ///< state speed to evaluate jacobians
97 ChMatrixRef mK, ///< result dQ/dx
98 ChMatrixRef mR, ///< result dQ/dv
99 ChMatrixRef mM ///< result dQ/da
100 ) = 0;
101
102 /// Access the jacobians (if any, i.e. if this is a stiff load)
GetJacobians()103 ChLoadJacobians* GetJacobians() { return this->jacobians; }
104
105 /// Create the jacobian loads if needed, and also
106 /// set the ChVariables referenced by the sparse KRM block.
107 virtual void CreateJacobianMatrices() = 0;
108
109 /// Update: this is called at least at each time step.
110 /// - It recomputes the generalized load Q vector(s)
111 /// - It recomputes the jacobian(s) K,R,M in case of stiff load
112 /// Q and jacobians assumed evaluated at the current state.
113 /// Jacobian structures are automatically allocated if needed.
114 virtual void Update(double time);
115
116 /// Report if this is load is stiff. If so, InjectKRMmatrices will provide
117 /// the jacobians of the load.
118 virtual bool IsStiff() = 0;
119
120 //
121 // Functions for interfacing to the state bookkeeping and solver
122 //
123
124 /// Adds the internal loads Q (pasted at global nodes offsets) into
125 /// a global vector R, multiplied by a scaling factor c, as
126 /// R += forces * c
127 virtual void LoadIntLoadResidual_F(ChVectorDynamic<>& R, const double c) = 0;
128
129 /// Tell to a system descriptor that there are item(s) of type
130 /// ChKblock in this object (for further passing it to a solver)
131 /// Basically does nothing, but inherited classes must specialize this.
132 virtual void InjectKRMmatrices(ChSystemDescriptor& mdescriptor);
133
134 /// Adds the current stiffness K and damping R and mass M matrices in encapsulated
135 /// ChKblock item(s), if any. The K, R, M matrices are added with scaling
136 /// values Kfactor, Rfactor, Mfactor.
137 virtual void KRMmatricesLoad(double Kfactor, double Rfactor, double Mfactor);
138 };
139
140 // -----------------------------------------------------------------------------
141
142 /// Class for a load acting on a single ChLoadable item, via ChLoader objects.
143 /// There are various ChLoader interfaces ready to use, that can be used
144 /// as 'building blocks'. These are especially important for creating loads
145 /// that are distributed on surfaces, lines, volumes, since some ChLoaders implement quadrature.
146 /// Create them as ChLoad< ChLoaderPressure > my_load(...); for example.
147
148 template <class Tloader>
149 class ChLoad : public ChLoadBase {
150 public:
151 Tloader loader;
152
ChLoad(std::shared_ptr<typename Tloader::type_loadable> mloadable)153 ChLoad(std::shared_ptr<typename Tloader::type_loadable> mloadable) : loader(mloadable) {}
154
~ChLoad()155 virtual ~ChLoad() {}
156
157 /// "Virtual" copy constructor (covariant return type).
Clone()158 virtual ChLoad* Clone() const override { return new ChLoad(*this); }
159
160 virtual int LoadGet_ndof_x() override;
161 virtual int LoadGet_ndof_w() override;
162 virtual void LoadGetStateBlock_x(ChState& mD) override;
163 virtual void LoadGetStateBlock_w(ChStateDelta& mD) override;
164 virtual void LoadStateIncrement(const ChState& x, const ChStateDelta& dw, ChState& x_new) override;
165 virtual int LoadGet_field_ncoords() override;
166
167 /// Compute Q, the generalized load.
168 /// Q is stored in the wrapped ChLoader.
169 /// Called automatically at each Update().
170 virtual void ComputeQ(ChState* state_x, ///< state position to evaluate Q
171 ChStateDelta* state_w ///< state speed to evaluate Q
172 ) override;
173
174 /// Compute jacobians (default fallback).
175 /// Uses a numerical differentiation for computing K, R, M jacobians, if stiff load.
176 /// If possible, override this with an analytical jacobian.
177 /// Compute the K=-dQ/dx, R=-dQ/dv , M=-dQ/da jacobians.
178 /// Called automatically at each Update().
179 virtual void ComputeJacobian(ChState* state_x, ///< state position to evaluate jacobians
180 ChStateDelta* state_w, ///< state speed to evaluate jacobians
181 ChMatrixRef mK, ///< result dQ/dx
182 ChMatrixRef mR, ///< result dQ/dv
183 ChMatrixRef mM ///< result dQ/da
184 ) override;
185
186 virtual void LoadIntLoadResidual_F(ChVectorDynamic<>& R, const double c) override;
187
188 /// Default: load is stiff if the loader is stiff. Override if needed.
IsStiff()189 virtual bool IsStiff() override { return loader.IsStiff(); }
190
191 /// Create the jacobian loads if needed, and also
192 /// set the ChVariables referenced by the sparse KRM block.
193 virtual void CreateJacobianMatrices() override;
194 };
195
196 // -----------------------------------------------------------------------------
197
198 /// Loads acting on a single ChLoadable item.
199 /// Differently form ChLoad, this does not use the ChLoader interface,
200 /// so one must inherit from this and implement ComputeQ() directly.
201 /// The ComputeQ() must write the generalized forces Q into the "load_Q" vector of this object.
202
203 class ChApi ChLoadCustom : public ChLoadBase {
204 public:
205 std::shared_ptr<ChLoadable> loadable;
206 ChVectorDynamic<> load_Q;
207
208 ChLoadCustom(std::shared_ptr<ChLoadable> mloadable);
209
~ChLoadCustom()210 virtual ~ChLoadCustom() {}
211
212 virtual int LoadGet_ndof_x() override;
213 virtual int LoadGet_ndof_w() override;
214 virtual void LoadGetStateBlock_x(ChState& mD) override;
215 virtual void LoadGetStateBlock_w(ChStateDelta& mD) override;
216 virtual void LoadStateIncrement(const ChState& x, const ChStateDelta& dw, ChState& x_new) override;
217 virtual int LoadGet_field_ncoords() override;
218
219 /// Compute jacobians (default fallback).
220 /// Uses a numerical differentiation for computing K, R, M jacobians, if stiff load.
221 /// If possible, override this with an analytical jacobian.
222 /// Compute the K=-dQ/dx, R=-dQ/dv , M=-dQ/da jacobians.
223 /// Called automatically at each Update().
224 virtual void ComputeJacobian(ChState* state_x, ///< state position to evaluate jacobians
225 ChStateDelta* state_w, ///< state speed to evaluate jacobians
226 ChMatrixRef mK, ///< result dQ/dx
227 ChMatrixRef mR, ///< result dQ/dv
228 ChMatrixRef mM ///< result dQ/da
229 ) override;
230
231 virtual void LoadIntLoadResidual_F(ChVectorDynamic<>& R, const double c) override;
232
233 /// Create the jacobian loads if needed, and also
234 /// set the ChVariables referenced by the sparse KRM block.
235 virtual void CreateJacobianMatrices() override;
236
237 /// Access the generalized load vector Q.
GetQ()238 virtual ChVectorDynamic<>& GetQ() { return load_Q; }
239 };
240
241 // -----------------------------------------------------------------------------
242
243 /// Loads acting on multiple ChLoadable items.
244 /// One must inherit from this and implement ComputeQ() directly. The ComputeQ() must
245 /// write the generalized forces Q into the "load_Q" vector of this object.
246 /// Given that multiple ChLoadable objects are referenced here, their sub-forces Q are
247 /// assumed appended in sequence in the "load_Q" vector, in the same order that has been
248 /// used in the std::vector "mloadables" for ChLoadCustomMultiple creation.
249 /// The same applies for the order of the sub-matrices of jacobians K,R etc.
250
251 class ChApi ChLoadCustomMultiple : public ChLoadBase {
252 public:
253 std::vector<std::shared_ptr<ChLoadable>> loadables;
254 ChVectorDynamic<> load_Q;
255
256 ChLoadCustomMultiple(std::vector<std::shared_ptr<ChLoadable>>& mloadables);
257 ChLoadCustomMultiple(std::shared_ptr<ChLoadable> mloadableA, std::shared_ptr<ChLoadable> mloadableB);
258 ChLoadCustomMultiple(std::shared_ptr<ChLoadable> mloadableA,
259 std::shared_ptr<ChLoadable> mloadableB,
260 std::shared_ptr<ChLoadable> mloadableC);
261
~ChLoadCustomMultiple()262 virtual ~ChLoadCustomMultiple() {}
263
264 virtual int LoadGet_ndof_x() override;
265 virtual int LoadGet_ndof_w() override;
266 virtual void LoadGetStateBlock_x(ChState& mD) override;
267 virtual void LoadGetStateBlock_w(ChStateDelta& mD) override;
268 virtual void LoadStateIncrement(const ChState& x, const ChStateDelta& dw, ChState& x_new) override;
269 virtual int LoadGet_field_ncoords() override;
270
271 /// Compute jacobians (default fallback).
272 /// Compute the K=-dQ/dx, R=-dQ/dv , M=-dQ/da jacobians.
273 /// Uses a numerical differentiation for computing K, R, M jacobians, if stiff load.
274 /// If possible, override this with an analytical jacobian.
275 /// NOTE: Given that multiple ChLoadable objects are referenced here, sub-matrices of mK,mR are
276 /// assumed pasted in i,j block-positions where i,j reflect the same order that has been
277 /// used in the std::vector "mloadables" at ChLoadCustomMultiple creation.
278 /// Called automatically at each Update().
279 virtual void ComputeJacobian(ChState* state_x, ///< state position to evaluate jacobians
280 ChStateDelta* state_w, ///< state speed to evaluate jacobians
281 ChMatrixRef mK, ///< result dQ/dx
282 ChMatrixRef mR, ///< result dQ/dv
283 ChMatrixRef mM ///< result dQ/da
284 ) override;
285
286 virtual void LoadIntLoadResidual_F(ChVectorDynamic<>& R, const double c) override;
287
288 /// Create the jacobian loads if needed, and also
289 /// set the ChVariables referenced by the sparse KRM block.
290 virtual void CreateJacobianMatrices() override;
291
292 /// Access the generalized load vector Q.
GetQ()293 virtual ChVectorDynamic<>& GetQ() { return load_Q; }
294 };
295
296 // =============================================================================
297 // IMPLEMENTATION OF ChLoad<Tloader> methods
298 // =============================================================================
299
300 template <class Tloader>
LoadGet_ndof_x()301 inline int ChLoad<Tloader>::LoadGet_ndof_x() {
302 return this->loader.GetLoadable()->LoadableGet_ndof_x();
303 }
304
305 template <class Tloader>
LoadGet_ndof_w()306 inline int ChLoad<Tloader>::LoadGet_ndof_w() {
307 return this->loader.GetLoadable()->LoadableGet_ndof_w();
308 }
309
310 template <class Tloader>
LoadGetStateBlock_x(ChState & mD)311 inline void ChLoad<Tloader>::LoadGetStateBlock_x(ChState& mD) {
312 this->loader.GetLoadable()->LoadableGetStateBlock_x(0, mD);
313 }
314
315 template <class Tloader>
LoadGetStateBlock_w(ChStateDelta & mD)316 inline void ChLoad<Tloader>::LoadGetStateBlock_w(ChStateDelta& mD) {
317 this->loader.GetLoadable()->LoadableGetStateBlock_w(0, mD);
318 }
319
320 template <class Tloader>
LoadStateIncrement(const ChState & x,const ChStateDelta & dw,ChState & x_new)321 inline void ChLoad<Tloader>::LoadStateIncrement(const ChState& x, const ChStateDelta& dw, ChState& x_new) {
322 this->loader.GetLoadable()->LoadableStateIncrement(0, x_new, x, 0, dw);
323 }
324
325 template <class Tloader>
LoadGet_field_ncoords()326 inline int ChLoad<Tloader>::LoadGet_field_ncoords() {
327 return this->loader.GetLoadable()->Get_field_ncoords();
328 }
329
330 template <class Tloader>
ComputeQ(ChState * state_x,ChStateDelta * state_w)331 inline void ChLoad<Tloader>::ComputeQ(ChState* state_x, ChStateDelta* state_w) {
332 this->loader.ComputeQ(state_x, state_w);
333 }
334
335 template <class Tloader>
ComputeJacobian(ChState * state_x,ChStateDelta * state_w,ChMatrixRef mK,ChMatrixRef mR,ChMatrixRef mM)336 inline void ChLoad<Tloader>::ComputeJacobian(ChState* state_x,
337 ChStateDelta* state_w,
338 ChMatrixRef mK,
339 ChMatrixRef mR,
340 ChMatrixRef mM) {
341 double Delta = 1e-8;
342
343 int mrows_w = this->LoadGet_ndof_w();
344 int mrows_x = this->LoadGet_ndof_x();
345
346 // compute Q at current speed & position, x_0, v_0
347 ChVectorDynamic<> Q0(mrows_w);
348 this->loader.ComputeQ(state_x, state_w); // Q0 = Q(x, v)
349 Q0 = this->loader.Q;
350
351 ChVectorDynamic<> Q1(mrows_w);
352 ChVectorDynamic<> Jcolumn(mrows_w);
353 ChState state_x_inc(mrows_x, nullptr);
354 ChStateDelta state_delta(mrows_w, nullptr);
355
356 // Compute K=-dQ(x,v)/dx by backward differentiation
357 state_delta.setZero(mrows_w, nullptr);
358
359 for (int i = 0; i < mrows_w; ++i) {
360 state_delta(i) += Delta;
361 this->LoadStateIncrement(*state_x, state_delta,
362 state_x_inc); // exponential, usually state_x_inc(i) = state_x(i) + Delta;
363 this->loader.ComputeQ(&state_x_inc, state_w); // Q1 = Q(x+Dx, v)
364 Q1 = this->loader.Q;
365 state_delta(i) -= Delta;
366
367 Jcolumn = (Q1 - Q0) * (-1.0 / Delta); // - sign because K=-dQ/dx
368 this->jacobians->K.block(0, i, mrows_w, 1) = Jcolumn;
369 }
370 // Compute R=-dQ(x,v)/dv by backward differentiation
371 for (int i = 0; i < mrows_w; ++i) {
372 (*state_w)(i) += Delta;
373 this->loader.ComputeQ(state_x, state_w); // Q1 = Q(x, v+Dv)
374 Q1 = this->loader.Q;
375 (*state_w)(i) -= Delta;
376
377 Jcolumn = (Q1 - Q0) * (-1.0 / Delta); // - sign because R=-dQ/dv
378 this->jacobians->R.block(0, i, mrows_w, 1) = Jcolumn;
379 }
380 }
381
382 template <class Tloader>
LoadIntLoadResidual_F(ChVectorDynamic<> & R,const double c)383 inline void ChLoad<Tloader>::LoadIntLoadResidual_F(ChVectorDynamic<>& R, const double c) {
384 unsigned int rowQ = 0;
385 for (int i = 0; i < this->loader.GetLoadable()->GetSubBlocks(); ++i) {
386 if (this->loader.GetLoadable()->IsSubBlockActive(i)) {
387 unsigned int moffset = this->loader.GetLoadable()->GetSubBlockOffset(i);
388 for (unsigned int row = 0; row < this->loader.GetLoadable()->GetSubBlockSize(i); ++row) {
389 R(row + moffset) += this->loader.Q(rowQ) * c;
390 ++rowQ;
391 }
392 }
393 }
394 }
395
396 template <class Tloader>
CreateJacobianMatrices()397 inline void ChLoad<Tloader>::CreateJacobianMatrices() {
398 if (!this->jacobians) {
399 // create jacobian structure
400 this->jacobians = new ChLoadJacobians;
401 // set variables forsparse KRM block
402 std::vector<ChVariables*> mvars;
403 loader.GetLoadable()->LoadableGetVariables(mvars);
404 this->jacobians->SetVariables(mvars);
405 }
406 }
407
408 } // end namespace chrono
409
410 #endif
411