Mathematical objects in Chrono

Table of Contents

• Linear algebra
• Function objects
• Quadrature

This documentation component focuses on Chrono mathematical functions and classes. These concepts are quite ubiquitous in the rest of the Chrono API and they are discussed also in the Linear algebra API documentation.

Linear algebra

Handling vectors and matrices is a recurring theme throughout the Chrono API. Chrono uses Eigen3 for representing all matrices (dense and sparse) and vectors.

Dense matrices in Chrono are templated by the scalar type and have row-major storage order. All Chrono matrix and vector types below are simply aliases to appropriate Eigen matrix types; see Linear algebra and the ChMatrix.h header.

The ChVector and ChVector2 classes, used to represent 3D vectors in space and 2D vectors in a plane, respectively, are not Eigen types nor derived from Eigen matrices.

Matrices are indexed starting from 0, with (row,column) indexing:

\[ \mathbf{A}=\left[ \begin{array}{cccc} a_{0,0} & a_{0,1} & a_{0,2} & ... \\ a_{1,0} & a_{1,1} & a_{1,2} & ... \\ a_{2,0} & ... & ... & ... \\ a_{n_{rows}-1,0} & ... & ... & a_{n_{rows}-1,n_{cols}-1} \end{array} \right] \]

There are many matrix and vector specializations and some of their basic features are outlined next.

Dynamic size matrices. Use ChMatrixDynamic to create a matrix with generic size, say 12 rows x 4 columns. ChMatrixDynamic is templated by the scalar type, with double as the default.

chrono::ChMatrixDynamic<double> A(12,4);

Fixed size matrices. Use ChMatrixNM to create matrices that do not need to be resized and whose size is known at compile-time.

chrono::ChMatrixNM<double,4,4> B;

From the Eigen documentation:

When should one use fixed sizes and when should one prefer dynamic sizes? The simple answer is: use fixed sizes for very small sizes where you can, and use dynamic sizes for larger sizes or where you have to. For small sizes, especially for sizes smaller than (roughly) 16, using fixed sizes is hugely beneficial to performance, as it allows Eigen to avoid dynamic memory allocation and to unroll loops.

The limitation of using fixed sizes, of course, is that this is only possible when you know the sizes at compile time. Also, for large enough sizes, say for sizes greater than (roughly) 32, the performance benefit of using fixed sizes becomes negligible. Worse, trying to create a very large matrix using fixed sizes inside a function could result in a stack overflow, since Eigen will try to allocate the array automatically as a local variable, and this is normally done on the stack. Finally, depending on circumstances, Eigen can also be more aggressive trying to vectorize (use SIMD instructions) when dynamic sizes are used.

3x3 fixed size matrices. Use ChMatrix33 to create 3x3 matrices, which are mostly used to represent rotation matrices and 3D inertia tensors. ChMatrix33 is templated by the scalar type (with double as the default). This matrix type is derived from a 3x3 fixed-size Eigen matrix with row-major storage and offers several dedicated constructors and methods for coordinate and rotation operations.

chrono::ChMatrix33<> R;

Dynamic size column vectors. Use ChVectorDynamic to create a column vector (one-column matrix) with a generic number of rows.

chrono::ChVectorDynamic<double> v(12);

Fixed size column vectors. Use ChVectorN to create a column vector with fixed length (known at compile time).

chrono::ChVectorN<double,6> w;

Row vectors. Use ChRowVectorDynamic and ChRowVectorN to create row vectors (one-row matrices) with dynamic size and fixed size, respectively.

In addition to the above types, specialized 3x4, 4x3, and 4x4 matrices used in multibody formalism are defined in ChMatrixMBD.h.

Basic operations with matrices

Consult the Eigen API for all matrix and vector arithmetic operations, block operations, and linear system solution.

The following code (see demo_CH_linalg.cpp) illustrates basic operations with matrices.

    std::cout << "\n=== Creation and assignment ===\n" << std::endl;
    {
        // Matrix allocated on the heap.
        chrono::ChMatrixDynamic<> Md1(5, 7);
        chrono::ChMatrixDynamic<> Md2(4, 4);
 
        // Matrix allocated on the stack.
        chrono::ChMatrixNM<double, 4, 4> Ms;
 
        // 3x3 matrices, mostly used for coordinate transformations.
        chrono::ChMatrix33<> ma;
 
        Ms.setConstant(0.1);  // Fill a matrix with an element
        std::cout << Ms << std::endl;
 
        Md1.setRandom();         // initialize with random numbers
        Md1.transposeInPlace();  // transpose the matrix in place
        Md1.resize(2, 2);        // resize
        Md1.setZero();           // set all elements to zero
 
        Md2.setIdentity();  // Create a diagonal matrix
        Md2 *= 3;
 
        Md2(0, 0) = 10;  // Use the () operator to reference single elements
        Md2(1, 2) = Md2(0, 0) + Md2(1, 1);
        Md2(2, 2) = 4;
 
        chrono::ChMatrixDynamic<> M1(Md2);   // Copy constructor
        chrono::ChMatrixDynamic<> M2(-Ms);   // The - unary operator returns a negated matrix
        chrono::ChMatrixDynamic<> M3(8, 4);  // 8x4 uninitialized matrix
        M3 = M1;                             // Copy (resize as needed)
        M3 = M1.transpose();                 // transposed copy.
    }
 
    std::cout << "\n=== Matrix operations ===\n" << std::endl;
    {
        chrono::ChMatrixNM<double, 2, 3> A;
        chrono::ChMatrixDynamic<double> B(2, 3);
        A.setRandom();
        B.setRandom();
        chrono::ChMatrixDynamic<> result(A + B);
        result = A + B;
        result = A;
        result += B;
 
        // Different ways to do subtraction..
        result = A - B;
        result = A;
        result -= B;
 
        // Multiplication with scalar
        result = A * 10;
        result = 10 * A;
        std::cout << result << std::endl;
 
        result = A;
        result *= 10;
        std::cout << result << std::endl;
 
        // Matrix multiplications
        chrono::ChMatrixDynamic<> C(3, 2);
        chrono::ChMatrixDynamic<> D(2, 3);
        chrono::ChMatrixDynamic<double> CD = C * D;
        chrono::ChMatrixDynamic<double> CD_t = D.transpose() * C.transpose();
 
        chrono::ChMatrix33<> R(chrono::Q_from_AngX(chrono::CH_C_PI / 3));
        std::cout << "rot * matrix\n" << R * C << std::endl;
        std::cout << "matrix * rot\n" << D * R << std::endl;
    }
 
    std::cout << "\n=== Chrono extensions to Eigen::MatrixBase ===\n" << std::endl;
    {
        chrono::ChMatrixNM<double, 2, 3> A;
        A.setRandom();
        std::cout << "random 2x3 matrix A:\n" << A << std::endl;
        A.fillDiagonal(10.1);
        std::cout << "fill diagonal with 10.1:\n" << A << std::endl;
        A.fill(2.1);
        std::cout << "fill entire matrix with 2.1:\n" << A << std::endl;
 
        chrono::ChMatrixNM<double, 2, 3> B(A);
        B(1, 2) += 0.01;
        std::cout << "matrix B = A with B(1,2) incremented by 0.01\n" << B << std::endl;
        std::cout << "|A-B| < 0.1?   " << A.equals(B, 0.1) << std::endl;
        std::cout << "|A-B| < 0.001? " << A.equals(B, 0.001) << std::endl;
 
        chrono::ChVectorDynamic<> v(3);
        v << 2, 3, 4;
        chrono::ChVectorDynamic<> w(3);
        w << 0.1, 0.2, 0.3;
        std::cout << "||v||_wrms, w = " << v.wrmsNorm(w) << std::endl;
        std::cout << "||v + v||_wrms, w = " << (v + v).wrmsNorm(w) << std::endl;
 
        std::cout << "v + v + 1: " << (v + v) + 1 << std::endl;
        std::cout << "1 + v: " << 1 + v << std::endl;
    }
 
    std::cout << "\n=== Matrix comparison tests ===\n" << std::endl;
    {
        chrono::ChMatrixNM<double, 2, 3> A;
        A.setRandom();
        chrono::ChMatrixDynamic<> B(A);
        if (A == B) {
            std::cout << "Matrices are exactly equal" << std::endl;
        }
 
        // Tolerance comparison
        B(1, 2) += 0.001;
        if (A.equals(B, 0.002)) {
            std::cout << "Matrices are equal within tol 0.002 \n";
        }
    }
 
    std::cout << "\n=== Pasting matrices and matrix-blocks ===\n" << std::endl;
    {
        chrono::ChMatrixDynamic<> A(4, 6);
        A.setRandom();
        std::cout << A << "\n\n";
 
        chrono::ChMatrixNM<double, 2, 3> B;
        B << 1, 2, 3, 4, 5, 6;
        std::cout << B << "\n\n";
 
        chrono::ChMatrixNM<double, 5, 7> C;
        C.setZero();
        C.block<2, 3>(1, 1) << 10, 20, 30, 40, 50, 60;
        std::cout << C << "\n\n";
 
        chrono::ChMatrixDynamic<> X = A;
        X.block<2, 3>(1, 2) = B;
        std::cout << X << "\n\n";
 
        X.block<2, 3>(1, 2) += C.block<2, 3>(1, 1);
        std::cout << X << "\n\n";
    }
 
    std::cout << "\n=== 3x3 matrix times vector ===\n" << std::endl;
    {
        chrono::ChVector<> v(1, 2, 3);
        chrono::ChMatrix33<> A;
        A.setRandom();
 
        chrono::ChMatrix33<> B(1);
        std::cout << B << std::endl;
 
        // Vector transformation, typical product [A] * v
        chrono::ChVector<> res1 = A * v;
 
        // Inverse vector transformation, [A]' * v
        chrono::ChVector<> res2 = A.transpose() * v;
 
        std::cout << A << std::endl;
        std::cout << res1 << std::endl;
        std::cout << res2 << std::endl;
    }
 
    std::cout << "\n=== Custom 3x4, 4x3, and 4x4 matrices ===\n" << std::endl;
    {
        chrono::ChMatrix34<> G = Eigen::Matrix<double, 3, 4>::Ones(3, 4);
        chrono::ChQuaternion<> q(1, 2, 3, 4);
        chrono::ChVector<> v = G * q;
        std::cout << v << std::endl;
 
        chrono::ChMatrix43<> Gt = 2 * G.transpose();
        std::cout << Gt << std::endl;
        q = Gt * v;
        std::cout << q << std::endl;
        q = G.transpose() * v;
        std::cout << q << std::endl;
 
        chrono::ChMatrix33<> rot(chrono::Q_from_AngX(chrono::CH_C_PI / 6));
        std::cout << rot << std::endl;
        std::cout << rot.transpose() << std::endl;
 
        chrono::ChMatrix34<> res1 = rot * G;
        chrono::ChMatrix34<> res2 = rot * Gt.transpose();
        chrono::ChMatrix43<> res3 = Gt * rot;
        chrono::ChMatrix43<> res4 = G.transpose() * rot;
 
        chrono::ChMatrix44<> A44;
        A44.setRandom();
        chrono::ChQuaternion<> q2 = A44 * q;
 
        chrono::ChStarMatrix44<> X(chrono::ChQuaternion<>(1, 2, 3, 4));
        std::cout << "4x4 star matrix X:\n" << X << std::endl;
        X.semiTranspose();
        std::cout << "Semi-transpose X:\n" << X << std::endl;
        X.semiNegate();
        std::cout << "Semi-negate X:\n" << X << std::endl;
    }
 
    std::cout << "\n=== Use of ChTransform ===\n" << std::endl;
    {
        chrono::ChVector<> vl(2, 3, 4);        // local point to transform
        chrono::ChVector<> t(5, 6, 7);         // translation of coord system
        chrono::ChQuaternion<> q(1, 3, 4, 5);  // rotation of coord system (quaternion)
        q.Normalize();                         // as unit quaternion, must be normalized
        chrono::ChMatrix33<> R;                // rotation of coord system (rotation matrix)
        R.Set_A_quaternion(q);                 // set from quaternion
        chrono::Coordsys csys(t, q);           // coordinate system representing translation + rotation
 
        // Perform the transformation: v = t + [R] * v'
        // NOTE: all the following ways will give the same result, so you can use them equivalently!
        auto va1 = chrono::ChTransform<>::TransformLocalToParent(vl, t, R);
        auto va2 = chrono::ChTransform<>::TransformLocalToParent(vl, t, q);
        auto va3 = csys.TransformLocalToParent(vl);
        auto va4 = t + R * vl;
        std::cout << va2.Equals(va1, 1e-6) << " " << va3.Equals(va1, 1e-6) << " " << va4.Equals(va1, 1e-6) << std::endl;
 
        // Inverse transformation
        auto vl1 = chrono::ChTransform<>::TransformParentToLocal(va1, t, q);
        auto vl2 = chrono::ChTransform<>::TransformParentToLocal(va1, t, R);
        auto vl3 = csys.TransformParentToLocal(va1);
        auto vl4 = R.transpose() * (va1 - t);
        std::cout << vl1.Equals(vl, 1e-6) << " " << vl2.Equals(vl, 1e-6) << " " << vl3.Equals(vl, 1e-6) << " "
                  << vl4.Equals(vl, 1e-6) << std::endl;
    }
 
    std::cout << "\n=== Linear systems ===\n" << std::endl;
    {
        chrono::ChMatrixNM<double, 3, 3> A;
        Eigen::Vector3d b;
        A << 1, 2, 3, 4, 5, 6, 7, 8, 10;
        b << 3, 3, 4;
        std::cout << "matrix A:\n" << A << std::endl;
        std::cout << "vector b:\n" << b << std::endl;
        Eigen::Vector3d x = A.colPivHouseholderQr().solve(b);
        std::cout << "solution:\n" << x << std::endl;
        std::cout << "Ax-b:\n" << A * x - b << std::endl;
    }

Function objects

These ChFunction objects are used in many places in Chrono::Engine, and are used to represent y=f(x) functions, for example when introducing prescribed displacements in a linear actuator.

These functions are scalar,

\[ x \in \mathbb{R} \rightarrow y \in \mathbb{R} \]

and there are a predefined number of them that are ready to use, such as sine, cosine, constant, etc. If the predefined ones are not enough, the user can implement his custom function by inheriting from the base ChFunction class.

See chrono::ChFunction for API details and a list of subclasses.

Example 1

ChFunction_Ramp f_ramp;
 
f_ramp.Set_ang(0.1);    // set angular coefficient;
f_ramp.Set_y0(0.4);     // set y value for x=0;
 
// Evaluate y=f(x) function at a given x value, using Get_y() :
double y    = f_ramp.Get_y(10);
// Evaluate derivative df(x)/dx at a given x value, using Get_y_dx() :
double ydx  = f_ramp.Get_y_dx(10);
 
GetLog() << "   ChFunction_Ramp at x=0: y=" << y << "  dy/dx=" << ydx << "\n\n";

Example 2

Save values of a sine ChFunction into a file.

ChFunction_Sine f_sine;
 
f_sine.Set_amp(2);      // set amplitude;
f_sine.Set_freq(1.5);   // set frequency;
 
ChStreamOutAsciiFile file_f_sine ("f_sine_out.dat");
 
// Evaluate y=f(x) function along 100 x points, and its derivatives, 
// and save to file (later it can be loaded, for example, in Matlab)
for (int i=0; i<100; i++)
{
    double x = (double)i/50.0;
    double y = f_sine.Get_y(x);
    double ydx = f_sine.Get_y_dx(x);
    double ydxdx = f_sine.Get_y_dxdx(x);
    file_f_sine << x << " " << y << " " << ydx << " " << ydxdx << "\n";
}   

Example 3

Define a custom function.

The following class will be used as an example of how you can create custom functions based on the ChFunction interface.

There is at least one mandatory member function to implement: Get_y().

Note that the base class implements a default computation of derivatives Get_ydx() and Get_ydxdx() by using a numerical differentiation, however if you know the analytical expression of derivatives, you can override the base Get_ydx() and Get_ydxdx() too, for higher precision.

// First, define a custom class inherited from ChFunction
 
class ChFunction_MyTest : public ChFunction
{
public:
    ChFunction* new_Duplicate() {return new ChFunction_MyTest;} 
    double Get_y(double x) {return cos(x);} // just for test: simple cosine
};
 
ChFunction_MyTest f_test;
 
ChStreamOutAsciiFile file_f_test ("f_test_out.dat");
 
// Evaluate y=f(x) function along 100 x points, and its derivatives, 
// and save to file (later it can be loaded, for example, in Matlab)
for (int i=0; i<100; i++)
{
    double x = (double)i/50.0;
    double y = f_test.Get_y(x);
    double ydx = f_test.Get_y_dx(x);
    double ydxdx = f_test.Get_y_dxdx(x);
    file_f_test << x << " " << y << " " << ydx << " " << ydxdx << "\n";
}

Quadrature

Quadrature is an operation that computes integrals as when computing areas and volumes.

The following code shows how to use the Gauss-Legendre quadrature to compute the integral of a function
\( \mathbb{R} \mapsto \mathbb{R} \) over a 1D interval, or \( f: \mathbb{R}^2 \mapsto \mathbb{R}\) over a 2D interval or \( f: \mathbb{R}^3 \mapsto \mathbb{R}\) over a 3D interval:

\[ F_{1D}=\int^a_b f(x) dx \]

\[ F_{2D}=\int^{a_y}_{b_y}\int^{a_x}_{b_x} f(x,y) dx dy \]

\[ F_{3D}=\int^{a_z}_{b_z}\int^{a_y}_{b_y}\int^{a_x}_{b_x} f(x,y,z) dx dy dz \]

If the function is polynomial of degree N and the quadrature is of order N, the result is exact, otherwise it is approximate (using large N improves quality but remember that this type of integration is often used where N in the range 1..10 suffices, otherwise other integration methods might be better).

For N less than 10, the quadrature uses precomputed coefficients for maximum performance.

// Define a y=f(x) function by inheriting ChIntegrable1D:
class MySine1d : public ChIntegrable1D<double>
{
public: 
    void Evaluate (double& result, const double x)  {
        result = sin(x);
    }
};
 
// Create an object from the function class
MySine1d mfx;
 
// Invoke 6th order Gauss-Legendre quadrature on 0..PI interval:
double qresult; 
ChQuadrature::Integrate1D<double>(qresult, mfx,  0, CH_C_PI,  6);
 
GetLog()<< "Quadrature 1d result:" << qresult << " (analytic solution: 2.0) \n";
 
 
// Other quadrature tests, this time in 2D
 
class MySine2d : public ChIntegrable2D<double>
{
public: 
    void Evaluate (double& result, const double x, const double y) { result = sin(x); }
};
 
MySine2d mfx2d;
ChQuadrature::Integrate2D<double>(qresult, mfx2d, 0, CH_C_PI, -1,1, 6);
GetLog()<< "Quadrature 2d result:" << qresult << " (analytic solution: 4.0) \n";

Note that thanks to templating, one can also integrate m-dimensional (vectorial, tensorial) functions \( \mathbf{f}: \mathbb{R}^n \mapsto \mathbb{R}^m \) , for example:

\[ \mathbf{F}=\int^{a_y}_{b_y}\int^{a_x}_{b_x} \mathbf{f}(x,y) dx dy \quad \mathbf{F} \in \mathbb{R}^2 \]

class MySine2dM : public ChIntegrable2D< ChMatrixNM<double,2,1> >
{
public: 
    void Evaluate (ChMatrixNM<double,2,1>& result, const double x, const double y) 
    { 
        result(0) = x*y;
        result(1) = 0.5*y*y;
    }
};
 
MySine2dM mfx2dM;
ChMatrixNM<double,2,1> resultM;
ChQuadrature::Integrate2D< ChMatrixNM<double,2,1> >(resultM, mfx2dM, 0, 1, 0,3, 6);
GetLog()<< "Quadrature 2d matrix result:" << resultM << " (analytic solution: 2.25, 4.5) \n";