I wrote this code as I am auditing a course "Physically-based and biomechanical modeling of natural phenomena" by Dr. Prusinkiewicz at the University of Calgary (Jan-April 2015).
//Code written by Jey for demonstrating the biomechanics of branching bending.
//The original description and implementation is in the "L + C modeling language (called as LPFG) paper"by P. Prusinkiewicz, R. Karwowski, and B. Lane (2007).
//Also please refer the paper "Integrating biomechanics into developmental plant models expressed using L−systems" by C. Jirasek, P. Prusinkiewicz, and B. Moulia (2000). They implemented the 3D version of the same algorithm using the CPFG to model growing cherries, grapes, and cycad. The authors also incorporated tropism, secondary growth and collision detection (based on the penalty method).
// I implemented the code using the list to understand and learn more about STL containers for the 2D version. You can also can implement using the vector container.
//The code will compile only in the recent Version of L-studio as I am using features available in L-studio (e.g., zooming, panning, color palette, reading parameters from the view file.)
//Below can see the animation of the implementation and the code for the same. The stiffness of the inter nodes is kept low in these models. If you keep make it high, it behaves like a rod.
#include <vve.h>
#include <util/parms.h>
#include <util/palette.h>
#include <algorithms/insert.h>
#include <geometry/geometry.h>
#include <vector>
#include <list>
#include <storage/storage.h>
#include <storage/graph.h>
#include <storage/vector.h>
#include <vector>
#include <iterator>
#include <math.h>
using geometry::Point3d;
using util::norm;
using algorithms::insert;
using std::cout;
using std::endl;
struct interNodeContent
{
double mass; //Mass of a node
double sigmaMass; //Combined mass to the right of the node
double torque; //Combined torque
double elasticity; //Elasticity at a node
double interNodeLength; //Length of the internode
Point3d leftNodePos; //Position of the left node
Point3d pos; //Position of the node
Point3d H; //Orientation of the internode
bool visited;
bool selected;
int nodeId;
interNodeContent():mass(0), sigmaMass(0), torque(0), elasticity(0), interNodeLength(0), leftNodePos(0,0,0), pos(0,0,0), H(0,0,0), visited(false), selected(false), nodeId(0)
{}
};
// Class defining your model
// If you change the name, don't forget to change the last line of the file too
class MyModel : public Model
{
public:
util::Palette palette;
double mechSteps;
double mass;
double elasticity;
double interNodeLength;
double backColor;
int color1, color2, color3, colorFail, background;
double lineWidth;
double error;
double relax;
int Phase; //Phase variable for backward or forward propgation.
double maxError;
int N;
double minTorque;
double maxTorque;
std::list<interNodeContent> interNodes; //Create a doublly-linked list called interNodes
void readParms()
{
util::Parms parms("view.v");
parms("SimulationParameters", "mechSteps", mechSteps);
parms("SimulationParameters", "Phase", Phase);
parms("SimulationParameters", "N", N); //number of internodes.
parms("PhysicalParameters", "mass", mass);
parms("PhysicalParameters", "elasticity", elasticity);
parms("PhysicalParameters", "interNodeLength", interNodeLength);
parms("PhysicalParameters", "relax", relax);
parms("PhysicalParameters", "maxError", maxError);
parms("PhysicalParameters", "maxTorque", maxTorque);
parms("PhysicalParameters", "minTorque", minTorque);
parms("View", "BackColor", backColor);
}
// Here, reread the files when they are modified
void modifiedFiles( const std::set<std::string>& filenames )
{
forall(const std::string& fn, filenames)
{
if(fn == "pal.map")
palette.reread();
else if(fn == "view.v")
readParms();
}
}
MyModel(QObject *parent) : Model(parent)
, palette("pal.map")
{
readParms();
// Registering the configuration files
registerFile("pal.map");
registerFile("view.v");
Point3d prevPos(0,0,0);
for(int i=0;i<N;i++) //Generate the internode contents
{
interNodeContent idc;
idc.nodeId = i;
idc.mass = mass;
idc.elasticity = elasticity;
idc.leftNodePos = prevPos;
idc.sigmaMass = 0;
idc.H.x() = 1; idc.H.y()=idc.H.z()=0;
idc.pos = prevPos + interNodeLength*idc.H;
prevPos=idc.pos;
interNodes.push_back(idc);
}
//for(std::list<interNodeContent>::iterator it = interNodes.begin() ; it != interNodes.end(); ++it)
{
//cout<<"Node ID "<<it->nodeId<<endl;
}
}
Point3d vectorRotate(Point3d a, double alpha) //Rotate a vector by alpha (in radians)
{
Point3d c;
c.x() = a.x()*cos(alpha) - a.y()*sin(alpha);
c.y() = a.x()*sin(alpha) + a.y()*cos(alpha);
c.z()=0;
return c;
}
void propogateLeft() //Calculate cumulative torque in the backward direction
{
std::list<interNodeContent>::reverse_iterator it = interNodes.rbegin();
it++;
std::list<interNodeContent>::reverse_iterator prevIt, nextIt;
while(it != interNodes.rend())
{
prevIt = it;
prevIt--;
it->sigmaMass = it->mass + prevIt->sigmaMass;
double crossProduct = ((prevIt->H)^((it->sigmaMass)*Point3d(0,-9.81,0))).z();
it->torque = prevIt->torque + crossProduct;
++it;
}
for(std::list<interNodeContent>::iterator it = interNodes.begin() ; it != interNodes.end(); ++it)
{
if(it->torque>maxTorque)
maxTorque = it->torque;
if(it->torque<minTorque)
minTorque = it->torque;
}
//cout<<"Max torque "<<maxTorque<<endl;
//cout<<"Min torque "<<minTorque<<endl;
}
int propogateRight() //Update node angles and internode position by iterating the list in the forward direction
{
error = 0;
std::list<interNodeContent>::iterator it = interNodes.begin();
std::list<interNodeContent>::iterator prevIt;
it++;
while(it != interNodes.end())
{
prevIt = it;
prevIt--;
Point3d newEqulibriumVector = vectorRotate(prevIt->H, it->torque/it->elasticity);
Point3d differenceVector = newEqulibriumVector - it->H;
error+=util::norm(differenceVector);
it->H = (it->H + relax*differenceVector).normalize();
it->pos = prevIt->pos + interNodeLength*it->H;
++it;
}
if(error<maxError)
return 1;
else
return 0;
}
void step()
{
for(int i = 0; i<mechSteps; i++) //Mechsteps control the frame rate.
{
if(Phase==1) //Start with Phase==1 and propagate left
propogateLeft();
if(propogateRight()==1) //If propogateRight returns 0, then the equilbrium is established. Make Phase = 0.
{
Phase = 0;
}
else
{
//Go to Phase 1 (Propogate Left). As you can see from the above code, the propogateLeft ()
//and propogateRight() goes in alternating manner until propgateRight() returns 0;
Phase = 1;
}
}
}
void preDraw()
{
util::Palette::Color c = palette.getColor(background);
glClearColor(background, background, background, 1);
//glClearColor(1.0, 1.0, 1.0, 1.0);
glLineWidth(lineWidth);
}
void draw(Viewer* viewer)
{
Point3d pmin, pmax;
interNodeContent p = interNodes.front();
pmax = pmin = p.pos;
for(std::list<interNodeContent>::iterator it = interNodes.begin() ; it != interNodes.end(); ++it)
{
Point3d& pos = it->pos;
if(pos.x() < pmin.x())
pmin.x() = pos.x();
if(pos.x() > pmax.x())
pmax.x() = pos.x();
if(pos.y() < pmin.y())
pmin.y() = pos.y();
if(pos.y() > pmax.y())
pmax.y() = pos.y();
}
viewer->setSceneBoundingBox(Vec(pmin), Vec(pmax));
for(std::list<interNodeContent>::iterator it = interNodes.begin() ; it != interNodes.end(); ++it)
{
color1= (int) 255*fabs(it->torque-minTorque)/(fabs(maxTorque-minTorque));
palette.useColor(color1);
drawCircle(it->pos);
}
}
void drawCircle(const Point3d& pos)
{
glBegin(GL_TRIANGLE_FAN);
glVertex3f(pos.x(),pos.y(),pos.z()+0.2);
for(float i=0;i<2*M_PI;i+=0.1){
glVertex3f(0.5*cos(i)+pos.x(),0.5*sin(i)+pos.y(),pos.z()+0.2);
}
glVertex3f(0.5*cos(0)+pos.x(),0.5*sin(0)+pos.y(),pos.z()+0.2);
glEnd();
}
};
// Define 'MyModel' as the model to be used
DEFINE_MODEL(MyModel);
DEFINE_VIEWER(Viewer);