Exercise 20.4.1

Exercise 20.4.1#

Q1: Test the examples#

  1. Update the upsteam and get folder mechsystem

  2. Modify Cmakelist:

① First, add cmake version, project name and standard

cmake_minimum_required(VERSION 3.20)
project(mechsystem)
set (CMAKE_CXX_STANDARD 20)

② Then, include the head file in the root ‘src’

include_directories(${CMAKE_CURRENT_SOURCE_DIR})

add_executable (test_mass_spring mass_spring.cpp ${CMAKE_SOURCE_DIR}/src)

③ In the end, use pybind to create python module

pybind11_add_module(mass_spring
    bind_mass_spring.cpp
)

target_include_directories(mass_spring PRIVATE 
    ${CMAKE_CURRENT_SOURCE_DIR}
)

  1. Build and compile in .venv enviroment. (Since we install python in .venv enviroment in case of extra conflicts)

// activate the enviroment
source .venv/bin/activate

// build file
mkdir build
cd build

// use venv Python to open CMake
cmake -DPython_EXECUTABLE=$(which python) ..

// Compile 
cmake ..
make -j

  1. We can test the exmaple files:

./test_mass_spring > test.txt

Test1: The result link

cd ~/ASC-ODE-team09/mechsystem

PYTHONPATH=../build/mechsystem python3 test_mass_spring.py >test2.txt

Test2 : The result link

  1. We can also test the 3D version according to mass_spring.ipynb

① Install pythreejs and activate

pip install pythreejs

jupyter nbextension enable --py widgetsnbextension
jupyter nbextension enable --py pythreejs

② Here we run directly in vscode. Pay attention: change the local repo

sys.path.append('../build/mechsystem')

In the example code, the spring constant is too large such that we can’t see an obviously ‘spring’ action.

If we change the spring constant with ‘1000’ and ‘500’, we can see:

Q2: Add Lagragian constraints#

  1. In order to add distance constraints to the MassSpring system, which means:

\[ | x_{i} - x_{j} | = L_{0} \]

A Lagrange multiplier 𝜆 enforces this constraint, leading to the Lagrange function：

\[ L(x,𝜆) = -U(x) + 𝜆g(x) \]

The system changes from ODE to DAE:

\[\begin{split} \begin{pmatrix} M & G^T \\ G & 0 \end{pmatrix} \begin{pmatrix} \ddot{x} \\ \lambda \end{pmatrix} = \begin{pmatrix} F_\text{ext} \\ 0 \end{pmatrix}, \end{split}\]

In the implementation, each constraint contributes:

  • a constraint equation \(g(x) = L-L_{0} = 0\)

  • a constraint force \(F_{i} ​= −λn, F_{j}​ = +λn\), n is unit direction vector

  1. Modify mass_spring.hpp and classify new function DistanceConstraint

class DistanceConstraint
{
public:
    Connector c1;
    Connector c2;
    double rest_length;
};

# Add new container
std::vector<DistanceConstraint> m_constraints;

# Add function
void addDistanceConstraint(const DistanceConstraint &dc)
{
    m_constraints.push_back(dc);
}

Mianly focus on MSS_Function:

for (size_t i = 0; i < n_constraints; i++)
    {
        auto& dc = mss.constraints()[i];
        double lambda = x(D * n_masses + i); 

        Vec<D> p1 = (dc.c1.type == Connector::FIX) ? mss.fixes()[dc.c1.nr].pos : xmat.row(dc.c1.nr);
        Vec<D> p2 = (dc.c2.type == Connector::FIX) ? mss.fixes()[dc.c2.nr].pos : xmat.row(dc.c2.nr);
        
        # Vector p2 -> p1
        Vec<D> d = p1 - p2; 
        double L = norm(d);
        if (L < 1e-12) continue;
        
        # Gradient of C with respect to p1
        Vec<D> dir = d / L; 

        # Constraint force
        if (dc.c1.type == Connector::MASS) fmat.row(dc.c1.nr) -= lambda * dir;
        if (dc.c2.type == Connector::MASS) fmat.row(dc.c2.nr) += lambda * dir;

        # Constraint Equation: g(x) = L - L0 = 0
        f(D * n_masses + i) = (L - dc.rest_length);
    }

Q3: Implement the exact derivative#

Since we are solving DAE, F(x) is non-linear. We need Jacobian derivative:

\[ J = \frac{∂F}{∂x} \]

The Jacobian has the block form

\[\begin{split} J = \begin{pmatrix} \frac{\partial F}{\partial x} & \frac{\partial F}{\partial \lambda} \\ \frac{\partial g}{\partial x} & 0 \end{pmatrix}, \end{split}\]

where

  • F(x) : physical forces (gravity + springs + constraint forces)

  • g(x)=0 : constraint equations

  • G = \nabla g(x) : gradient of the constraints

  • \(\lambda\): Lagrange multipliers (constraint forces)

  1. First, do the derivative of the spring forces. Taking the derivative yields the geometric stiffness matrix:

\[ K = k\, nn^T + \frac{k(L-L_0)}{L}(I - nn^T). \]
double ninj = n(i)*n(j);
double Kij = k * ninj + force_over_L * ((i==j?1.0:0.0) - ninj);

# Assemble into Jacobian df
df(c1.nr*D + i, c1.nr*D + j) -= Kij;
df(c2.nr*D + i, c2.nr*D + j) -= Kij;
df(c1.nr*D + i, c2.nr*D + j) += Kij;
df(c2.nr*D + i, c1.nr*D + j) += Kij;

  1. Then, do the derivative of the constraint forces:

# Geometric stiffness due to constraint tension (lambda)

double lambda_over_L = lambda / L;


# Hessian of the constraint: (lambda/L) * (I - n*n^T)

if (dc.c1.type == Connector::MASS) df(dc.c1.nr*D + i, dc.c1.nr*D + j) -= Hij;
if (dc.c2.type == Connector::MASS) df(dc.c2.nr*D + i, dc.c2.nr*D + j) -= Hij;
if (dc.c1.type == Connector::MASS && dc.c2.type == Connector::MASS) 

df(dc.c1.nr*D + i, dc.c2.nr*D + j) += Hij;
df(dc.c2.nr*D + i, dc.c1.nr*D + j) += Hij;

  1. Finally, coupling between x and \(\lambda\)

df(dc.c1.nr*D + i, idx_lambda) -= n(i);  // dF/dlambda
df(idx_lambda, dc.c1.nr*D + i) += n(i);  // dC/dx

3D Version#

For the last question, we successfully simulate 4 different models, the codes and result is in next page.

  • Extend the double-pendulum to a chain

  • Build more complex structures using beams (for example a crane), and simulate vibration

  • Replace the springs of the double-pendulum by distance-constraints

  • Build a spinning top (German: Kreisel).

Contents