test/fem-springs.cpp

#include <iostream>
#include <fstream>
#include <cmath>

/* TYPE AND FUNCTION DEFINITIONS */
struct Vertex;
struct Element;
void doCompute();

struct Vertex {
   Vertex(int    _id,
          double _xy[2],
         double _m) {
      id = _id;
      xy[0] = _xy[0]; xy[1] = _xy[1];
      u0[0] = 0.0; u0[1] = 0.0;
      u1[0] = 0.0; u1[1] = 0.0;
      f[0] = 0.0; f[1] = 0.0;
      m = _m; }
   
   int id;
   double xy[2];
   double u0[2];
   double u1[2];
   double f[2];
   double m;
};

struct Element {
   int id;
   Vertex*  vertices[2];
   Element* adjacents[3];
};

/* GLOBAL VARIABLES */
double dt = 0.0001;
double t = 0.0;
double end_time = 0.5;
double k = 1000.0;
double d = 5.0;
Vertex va(0, (double[]){0,0}, 0.1);
Vertex vb(1, (double[]){0,-5}, 0.1);

int main() {
   // Initialise file to write to
   std::ofstream outfile;
   outfile.open("fem-springs-output.txt");
   
   // Set initial conditions
   vb.u0[0] = 0.0;
   vb.u0[1] = -1.0;
   vb.u1[0] = 0.0;
   vb.u1[1] = -1.0;
   
   std::cout.setf(std::ios::fixed,std::ios::floatfield);
   std::cout.precision(8);
   outfile.setf(std::ios::fixed,std::ios::floatfield);
   outfile.precision(8);
   
   std::cout << "T\t\tX\t\tY" << std::endl;
   outfile   << "T\t\tX\t\tY" << std::endl;
   while(t < end_time) {
      doCompute();
      std::cout << t << "\t" << vb.xy[0] + vb.u1[0] << "\t" << vb.xy[1] + vb.u1[1] << std::endl;
      outfile   << t << "\t" << vb.xy[0] + vb.u1[0] << "\t" << vb.xy[1] + vb.u1[1] << std::endl;
      t+=dt;
   }
   
   // Close output file
   outfile.close();
   
   return 0;
}

void doCompute() {
   // Compute forces
   double A[2], a[2], b[2];
   for(int i=0; i<2; i++) {
      A[i] = vb.xy[i] - va.xy[i];
      a[i] = A[i] + vb.u1[i] - va.u1[i];
      b[i] = A[i] + vb.u0[i] - va.u0[i];
   }
   
   double modA = sqrt((A[0] * A[0]) + (A[1] * A[1]));
   double moda = sqrt((a[0] * a[0]) + (a[1] * a[1]));
   double modb = sqrt((b[0] * b[0]) + (b[1] * b[1]));
   
   for(int i=0; i<2; i++) {
      vb.f[i] = (k * (a[i] - (a[i]*(modA/moda)))) + (d * ((a[i] - (a[i]*(modb/moda)))/dt));
   }
   
   // Compute displacements
   double u2[2];
   for(int i=0; i<2; i++) {
      u2[i] = (((vb.m*(2*vb.u1[i]-vb.u0[i]))/(dt*dt))-vb.f[i])*((dt*dt)/vb.m);
   }
   
   // Update displacements arrays
   for(int i=0; i<2; i++) {
      vb.u0[i] = vb.u1[i];
      vb.u1[i] = u2[i];
   }
}

---