masonium.com

CS 176 Project 4-5: Discrete Elastic Rods

Introduction

Much of computer graphics is dedicated to the accurate simulation of physical phenomena. Previous attempts at simulating elastic materials did not accurately model certain aspects. Using the techniques explored in the Discrete Elastic Rods paper, we can more accurately reproduce these phenomena

My program uses the bending and twisting energies to find a quasi-static equilibrium under inextensibility constraints.

Running the Program

You can download the source here. Compiling the program require OpenGL, GLUT, Boost, Boost.ProgramOptions, BLAS, and LAPACK. A Makefile is included, so building the program should be as simple as `make`.

You can download a Linux binary here.

The program options are as follows:

Options:
  -i [ --input-file ] arg          input file (required)
  -s [ --step-size ] arg (=0.0001) step size to use for integration
  -e [ --euler ]                   use an euler solver, rather than
                                   the default runge-kutta solver
  -h [ --help ]                    display help

Elastic Rod Options:
  -s [ --stiffness ] arg (=0.1)     stiffness constant \\alpha
  -t [ --twist-modulus ] arg (=0.1) twisting modulus \\beta
  -i [ --initial-twist ] arg (=0)   initial amount of total 
                                    twist in the rod between 
                                    the two endpoints

The simulation starts as soon as the program does.

Implementation

File Format

25
1.00000 0.00000 0.00000 0
0.96593 0.25882 0.00000 1
0.86603 0.50000 0.00000 1
0.70711 0.70711 0.00000 1
0.50000 0.86603 0.00000 1
0.25882 0.96593 0.00000 1
0.000000 1.00000 0.00000 1
-0.25882 0.96593 0.00000 1
...

The first number is just the number of points. Each subsequent line represents a point. The first three numbers represent the (x, y, z) coordinates, and the last number indicates whether or not the point is free to move or is fixed. (A '1' represents a free point).

Inextensibility Constraint

The inextensibility force direction is computed as the negative gradient of the stretch energy. Instead of using the magnitude of this force, it is scaled to be the a multiple of the magnitude of the total bending and twisting force.

Results

Bending Only

For these first examples, the left side is the initial phase, and the right side is final output.

In this first example, the endpoints of the sin wave are fixed.

In this example, the first and last endpoints are the same point in space and are both fixed.

Bending and Twisting

(Left): We add a small initial twist angle to the endpoints and view the ring from the side. The ring actually begins to rotate after its rigid shape stabilizes. (Middle, Right): When we add a larger initial twist, the ring folds in on itself.

(Left): The parabola has its two endpoints fixed. (Middle): The parabola begins to curl up, since there is no gravity. (Right): Plectonemes begin to form in the rod.

Discussion

Since the program does not implement collision detection, twisting rods rarely stabilize. Rather, the folds tend to twist periodically.