git.sesse.net Git - sproing/blob - main.cpp

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <stdlib.h>
   4 #include <unistd.h>
   5 #include <sys/ioctl.h>
   6 #include <vector>
   7 #include <algorithm>
   8 #include "glwindow.h"
   9
  10 struct point {
  11         double x, y;
  12         double px, py;
  13         double ax, ay;
  14 };
  15 struct measurement {
  16         int p1, p2;
  17         double d;
  18         double stress;
  19         double alpha;
  20 };
  21 struct angle {
  22         int p1, p2, p3;
  23         double rad;
  24         double stress;
  25         double alpha;
  26 };
  27 std::vector<point> points;
  28 std::vector<measurement> measurements;
  29 std::vector<angle> angles;
  30
  31 void reload()
  32 {
  33         FILE *file = fopen("spec.txt", "r");
  34         unsigned num_points;
  35         fscanf(file, "%u", &num_points);
  36
  37         if (num_points < points.size()) {
  38                 points.erase(points.begin() + num_points, points.end());
  39         } else if (num_points > points.size()) {
  40                 for (int i = points.size(); i < num_points; ++i) {
  41                         point p;
  42                         if (i == 0) {
  43                                 p.x = p.y = 0.0;
  44                         } else {
  45                                 p.x = rand() / (RAND_MAX + 1.0);
  46                                 p.y = rand() / (RAND_MAX + 1.0);
  47                         }
  48                         if (i == 1) {
  49                                 p.y = 0.0;
  50                         }
  51                         p.px = p.x;
  52                         p.py = p.y;
  53                         points.push_back(p);
  54                 }
  55         }
  56
  57         measurements.erase(measurements.begin(), measurements.end());
  58         angles.erase(angles.begin(), angles.end());
  59
  60         for ( ;; ) {
  61                 char buf[256];
  62                 if (fscanf(file, "%s", &buf) != 1)
  63                         break;
  64
  65                 if (buf[0] == 'd') {
  66                         measurement l;
  67                         if (fscanf(file, "%u %u %lf %lf", &l.p1, &l.p2, &l.d, &l.alpha) != 4)
  68                                 break;
  69                         measurements.push_back(l);
  70                 }
  71                 if (buf[0] == 'a') {
  72                         angle a;
  73                         if (fscanf(file, "%u %u %u %lf %lf", &a.p1, &a.p2, &a.p3, &a.rad, &a.alpha) != 5)
  74                                 break;
  75                         a.rad *= M_PI / 180.0;
  76                         angles.push_back(a);
  77                 }
  78         }
  79
  80         fclose(file);
  81 }
  82
  83 void recalc()
  84 {
  85         double dt = 1.0/60.0;
  86         double k = 0.01;
  87         double k_angle = 0.02;
  88         double f = 0.02 * dt;
  89
  90         for (unsigned i = 0; i < points.size(); ++i) {
  91                 points[i].ax = points[i].ay = 0.0;
  92         }
  93
  94         for (unsigned i = 0; i < measurements.size(); ++i) {
  95                 double dx = points[measurements[i].p2].x - points[measurements[i].p1].x;
  96                 double dy = points[measurements[i].p2].y - points[measurements[i].p1].y;
  97                 double d = sqrt(dx*dx + dy*dy);
  98                 measurements[i].stress = fabs(d - measurements[i].d);
  99
 100                 points[measurements[i].p1].ax += (d - measurements[i].d) * dx/d * k;
 101                 points[measurements[i].p1].ay += (d - measurements[i].d) * dy/d * k;
 102                 points[measurements[i].p2].ax -= (d - measurements[i].d) * dx/d * k;
 103                 points[measurements[i].p2].ay -= (d - measurements[i].d) * dy/d * k;
 104         }
 105
 106         for (unsigned i = 0; i < angles.size(); ++i) {
 107                 double cx = points[angles[i].p1].x - points[angles[i].p2].x;
 108                 double cy = points[angles[i].p1].y - points[angles[i].p2].y;
 109                 double c = sqrt(cx*cx + cy*cy);
 110
 111                 double bx = points[angles[i].p3].x - points[angles[i].p2].x;
 112                 double by = points[angles[i].p3].y - points[angles[i].p2].y;
 113                 double b = sqrt(bx*bx + by*by);
 114
 115                 double alpha = acos((bx*cx + by*cy) / (b*c));
 116                 angles[i].stress = fabs(alpha - angles[i].rad);
 117
 118                 double tp3x = -by / b;
 119                 double tp3y = bx / b;
 120
 121                 double tp1x = cy / c;
 122                 double tp1y = -cx / c;
 123
 124                 points[angles[i].p1].ax -= (alpha - angles[i].rad) * tp1x * k_angle;
 125                 points[angles[i].p1].ay -= (alpha - angles[i].rad) * tp1y * k_angle;
 126                 points[angles[i].p3].ax -= (alpha - angles[i].rad) * tp3x * k_angle;
 127                 points[angles[i].p3].ay -= (alpha - angles[i].rad) * tp3y * k_angle;
 128         }
 129
 130         // constraints
 131         if (points.size() >= 1) {
 132                 points[0].ax = points[0].ay = 0.0;
 133         }
 134         if (points.size() >= 2) {
 135                 points[1].ay = 0.0;
 136         }
 137
 138         // Verlet integration
 139         for (unsigned i = 0; i < points.size(); ++i) {
 140                 double nx = (2.0 - f) * points[i].x - (1.0 - f) * points[i].px + points[i].ax * dt * dt;
 141                 double ny = (2.0 - f) * points[i].y - (1.0 - f) * points[i].py + points[i].ay * dt * dt;
 142                 points[i].px = points[i].x;
 143                 points[i].py = points[i].y;
 144                 points[i].x = nx;
 145                 points[i].y = ny;
 146         }
 147
 148         if (points.size() >= 2) {
 149                 if (points[1].x < 0.0) {
 150                         points[1].x = -points[1].x;
 151                         points[1].px = -points[1].px;
 152                         points[1].ax = -points[1].ax;
 153                 }
 154         }
 155 }
 156
 157 void redraw(GLWindow &win)
 158 {
 159         glClear(GL_COLOR_BUFFER_BIT);
 160
 161         glEnable(GL_BLEND);
 162         glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
 163         //glEnable(GL_LINE_SMOOTH);
 164         //glLineWidth(5.0);
 165
 166         glMatrixMode(GL_PROJECTION);
 167         glLoadIdentity();
 168         glOrtho(-1.0, 1.0, -1.0, 1.0, 0.0, 1.0);
 169
 170         glMatrixMode(GL_MODELVIEW);
 171         glLoadIdentity();
 172         glScalef(0.1, 0.1, 0.1);
 173
 174         glPointSize(5.0);
 175         glColor4f(1.0, 1.0, 1.0, 1.0);
 176
 177         glBegin(GL_POINTS);
 178         for (int i = 0; i < points.size(); ++i) {
 179                 glVertex2f(points[i].x, points[i].y);
 180         }
 181         glEnd();
 182
 183         glBegin(GL_LINES);
 184         for (int i = 0; i < measurements.size(); ++i) {
 185                 double sfac = std::min(5.0 * (fabs(measurements[i].stress) / measurements[i].d), 1.0);
 186                 double r, g, b;
 187                 if (sfac < 0.5) {
 188                         r = 2.0 * sfac;
 189                         g = 1.0;
 190                         b = 0.0;
 191                 } else {
 192                         r = 1.0;
 193                         g = 1.0 - 2.0 * (sfac - 0.5);
 194                         b = 0.0;
 195                 }
 196
 197                 glColor4f(r, g, b, measurements[i].alpha);
 198                 glVertex2f(points[measurements[i].p1].x, points[measurements[i].p1].y);
 199                 glVertex2f(points[measurements[i].p2].x, points[measurements[i].p2].y);
 200         }
 201         glEnd();
 202
 203         win.flip();
 204 }
 205
 206 void randomize(double Z)
 207 {
 208         for (int i = 2; i < points.size(); ++i) {
 209                 points[i].x += Z * (rand() / (RAND_MAX + 1.0) - .5);
 210                 points[i].y += Z * (rand() / (RAND_MAX + 1.0) - .5);
 211         }
 212 }
 213
 214 void destress()
 215 {
 216         double max_stress = 0.0;
 217         int p1, p2;
 218
 219         for (int i = 0; i < measurements.size(); ++i) {
 220                 if (measurements[i].p1 <= 1 ||
 221                     measurements[i].p2 <= 1)
 222                         continue;
 223
 224                 if (measurements[i].stress > max_stress) {
 225                         max_stress = measurements[i].stress;
 226                         p1 = measurements[i].p1;
 227                         p2 = measurements[i].p2;
 228                 }
 229         }
 230
 231         std::swap(points[p1], points[p2]);
 232 }
 233
 234 int main(void)
 235 {
 236         GLWindow glw("foo", 640, 480, 32, false, 16, -1);
 237         int one = 1;
 238
 239         srand((time_t)time(NULL));
 240
 241         ioctl(0, FIONBIO, &one);
 242         reload();
 243
 244         for ( ;; ) {
 245                 char ch;
 246                 while (read(0, &ch, 1) == 1) {
 247                         if (ch == 'r') {
 248                                 reload();
 249                         }
 250                         if (ch == 'z') {
 251                                 randomize(0.005);
 252                         }
 253                         if (ch == 'Z') {
 254                                 randomize(0.1);
 255                         }
 256                         if (ch == 's') {
 257                                 destress();
 258                         }
 259                 }
 260                 recalc();
 261                 recalc();
 262                 recalc();
 263                 recalc();
 264                 recalc();
 265                 redraw(glw);
 266         }
 267 }