strain_GL.cpp

//
// strain_GL.cpp
//    input strain, display GL
//    (c) MIT CBA Neil Gershenfeld 7/1/20
//
// includes
//
#include <iostream>
#include <cstring>
#include <GL/glew.h>
#include <GLFW/glfw3.h>
using namespace std;
//
// defines
//
#define CPP (8*3*6) // render coordinates per particle
   // 8 triangles * 3 verts * (3 coords + 3 colors)
#define VPP (8*3) // render vertices per particle
   // 8 triangles * 3 verts
//
// global variables
//
class vars {
   public:
      float rx; // view rotation
      float ry;
      float rz;
      float dr = 0.01;
      float tx = 0; // view translation
      float ty = 0;
      float tz = 0;
      float dt = 0.01;
      float sxyz; // view scale
      float scale; // strain scale multiplier
      float perspective; // view perspective
      int mx,my; // mouse position
      int mb = -1; // mouse state
      int np; // number of particles
      float xmin,xmax,ymin,ymax,zmin,zmax; // data limits
      float *px,*py,*pz; // particle positions
      char *type; // particle type
      float *s; // strain
      float d; // particle pitch
      float size; // particle size
      float* data; // GL drawing buffer
      GLuint program; // GL program
      int width = 1000; // GL window width
      int height = 1000; // GL window height
   } v;
//
// GL class
//
class GL {
   public:
      //
      // variables
      //
      GLFWwindow* window;
      GLuint vertex;
      GLuint fragment;
      GLuint buffer;
      GLint location;
      //
      // keyboard callback (escape to quit)
      //
      static void keyboard(GLFWwindow* window,int key,int scancode,int action,int mods) {
         GLint location;
         if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS) {
            glfwDestroyWindow(window);
            glfwTerminate();
            cerr << endl;
            exit(0);
            }
         else if (key == 'Q' && action == GLFW_PRESS) {
            glfwDestroyWindow(window);
            glfwTerminate();
            cerr << endl;
            exit(0);
            }
         else if (key == '-' && action == GLFW_PRESS) {
            v.sxyz /= 1.1;
            location = glGetUniformLocation(v.program,"sxyz");
            glUniform1f(location,v.sxyz);
            }
         else if (key == '=' && action == GLFW_PRESS) {
            v.sxyz *= 1.1;
            location = glGetUniformLocation(v.program,"sxyz");
            glUniform1f(location,v.sxyz);
            }
         }
      //
      // scroll callback
      //
      static void scroll(GLFWwindow* window,double x,double y) {
         GLint location;
         if (y > 0) {
            v.sxyz *= 1.1;
            location = glGetUniformLocation(v.program,"sxyz");
            glUniform1f(location,v.sxyz);
            }
         else {
            v.sxyz /= 1.1;
            location = glGetUniformLocation(v.program,"sxyz");
            glUniform1f(location,v.sxyz);
            }
         cerr << '\r' << "   scale " << v.sxyz << "   " << flush;
         }
      //
      // mouse cursor callback (left down to rotate)
      //
      static void cursor(GLFWwindow* window,double x,double y) {
         GLint location;
         if (v.mb == GLFW_MOUSE_BUTTON_LEFT) {
            if (v.mx < x)
               v.ry -= v.dr;
            else if (v.mx > x)
               v.ry += v.dr;
            v.mx = x;
            if (v.my < y)
               v.rx -= v.dr;
            else if (v.my > y)
               v.rx += v.dr;
            v.my = y;
            cerr << '\r' << "   rotate " << v.rx << ' ' << v.ry << "   " << flush;
            location = glGetUniformLocation(v.program,"rx");
            glUniform1f(location,v.rx);
            location = glGetUniformLocation(v.program,"ry");
            glUniform1f(location,v.ry);
            }
         }
      //
      // mouse button callback
      //
      static void button(GLFWwindow* window,int button,int action,int mods) {
         double x,y;
         if ((button == GLFW_MOUSE_BUTTON_LEFT)
               && (action == GLFW_PRESS)) {
            glfwGetCursorPos(window,&x,&y);
            v.mb = GLFW_MOUSE_BUTTON_LEFT;
            v.mx = x;
            v.my = y;
            }
         else if ((button == GLFW_MOUSE_BUTTON_LEFT)
               && (action == GLFW_RELEASE))
            v.mb = -1;
         }
      //
      // update and draw buffer
      //
      void draw(void (*update)()) {
         glfwPollEvents();
         update();
         glBufferSubData(GL_ARRAY_BUFFER,0,CPP*v.np*sizeof(float),v.data);
         glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);  
         glDrawArrays(GL_TRIANGLES,0,VPP*v.np);
         glfwSwapBuffers(window);
         }
      //
      // constructor
      //
      GL(int width,int height,char *title) {
         //
         // init GLFW
         //
         if (!glfwInit())
            cerr << "   could not start GLFW" << endl;
         glfwWindowHint(GLFW_SAMPLES,4);
         window = glfwCreateWindow(width,height,title,NULL,NULL);
         if (!window)
            cerr << "   could not open window" << endl;
         glfwMakeContextCurrent(window);
         glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
         glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
         glfwWindowHint(GLFW_OPENGL_PROFILE,GLFW_OPENGL_CORE_PROFILE);
         //
         // init GLEW
         //
         glewExperimental = GL_TRUE;
         glewInit();
         const GLubyte* renderer = glGetString(GL_RENDERER);
         const GLubyte* version = glGetString(GL_VERSION);
         cerr << "   renderer: " << renderer << endl;
         cerr << "   GL version " << version << endl;
         //
         // vertex shader
         //
         const char* vertex_code = R"(
            attribute vec3 position;
            attribute vec3 color;
            varying vec3 fcolor;
            uniform float rx,ry,rz,tx,ty,tz,sxyz,perspective;
            mat4 rotatex(float angle) {
               return mat4(1.0,0.0,0.0,0.0,
               0.0,cos(angle),-sin(angle),0.0,
               0.0,sin(angle),cos(angle),0.0,
               0.0,0.0,0.0,1.0);
               }
            mat4 rotatey(float angle) {
               return mat4(cos(angle),0.0,sin(angle),0.0,
               0.0,1.0,0.0,0.0,
               -sin(angle),0.0,cos(angle),0.0,
               0.0,0.0,0.0,1.0);
               }
            mat4 rotatez(float angle) {
               return mat4(cos(angle),-sin(angle),0.0,0.0,
               sin(angle),cos(angle),0.0,0.0,
               0.0,0.0,1.0,0.0,
               0.0,0.0,0.0,1.0);
               }
            mat4 translate(float dx,float dy,float dz) {
               return mat4(1.0,0.0,0.0,dx,
               0.0,1.0,0.0,dy,
               0.0,0.0,1.0,dz,
               0.0,0.0,0.0,1.0);
               }
            mat4 scale(float s) {
               return mat4(s,0,0.0,0.0,
               0,s,0.0,0.0,
               0.0,0.0,s,0.0,
               0.0,0.0,perspective,1.0);
               }
            void main() {
               fcolor = color;
               vec4 vertex = vec4(position,1.0);
               gl_Position = vertex*rotatex(rx)
                  *rotatey(ry)*rotatez(rz)
                  *translate(tx,ty,-.9*tz)*scale(sxyz);
               }
            )";
         //
         // fragment shader
         //
         const char* fragment_code = R"(
            varying vec3 fcolor;
            float depth;
            void main() {
               depth = 1.0-1.0*gl_FragCoord.z;
               gl_FragColor = vec4(depth*fcolor[0],depth*fcolor[1],depth*fcolor[2],1.0);
               }
            )";
         //
         // compile and link shaders
         //
         v.program = glCreateProgram();
         GLuint vertex = glCreateShader(GL_VERTEX_SHADER);
         GLuint fragment = glCreateShader(GL_FRAGMENT_SHADER);
         glShaderSource(vertex,1,&vertex_code,NULL);
         glShaderSource(fragment,1,&fragment_code,NULL);
         glCompileShader(vertex);
         glCompileShader(fragment);
         glAttachShader(v.program,vertex);
         glAttachShader(v.program,fragment);
         glLinkProgram(v.program);
         glDetachShader(v.program,vertex);
         glDetachShader(v.program,fragment);
         glUseProgram(v.program);
         //
         // set up buffers
         //
         glGenBuffers(1,&buffer);
         glBindBuffer(GL_ARRAY_BUFFER,buffer);
         glBufferData(GL_ARRAY_BUFFER,CPP*v.np*sizeof(float),
            v.data,GL_DYNAMIC_DRAW);
         location = glGetAttribLocation(v.program,"position");
         glEnableVertexAttribArray(location);
         glVertexAttribPointer(location,3,GL_FLOAT,GL_FALSE,
            6*sizeof(float),0);
         location = glGetAttribLocation(v.program,"color");
         glEnableVertexAttribArray(location);
         glVertexAttribPointer(location,3,GL_FLOAT,GL_FALSE,
            6*sizeof(float),(GLvoid*)(3*sizeof(float)));
         //
         // set up drawing
         //
         glEnable(GL_DEPTH_TEST);
         glDepthMask(GL_TRUE);
         glDepthFunc(GL_LESS);
         glDepthRange(0.0,1.0);
         glEnable(GL_BLEND);
         glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
         glEnable(GL_MULTISAMPLE);
         //
         // set up view
         //
         location = glGetUniformLocation(v.program,"rx");
         glUniform1f(location,v.rx);
         location = glGetUniformLocation(v.program,"ry");
         glUniform1f(location,v.ry);
         location = glGetUniformLocation(v.program,"rz");
         glUniform1f(location,v.rz);
         location = glGetUniformLocation(v.program,"tx");
         glUniform1f(location,v.tx);
         location = glGetUniformLocation(v.program,"ty");
         glUniform1f(location,v.ty);
         location = glGetUniformLocation(v.program,"tz");
         glUniform1f(location,v.tz);
         location = glGetUniformLocation(v.program,"sxyz");
         glUniform1f(location,v.sxyz);
         location = glGetUniformLocation(v.program,"perspective");
         glUniform1f(location,v.perspective);
         //
         // set up callbacks
         //
         glfwSetKeyCallback(window,keyboard);
         glfwSetScrollCallback(window,scroll);
         glfwSetMouseButtonCallback(window,button);
         glfwSetCursorPosCallback(window,cursor);
         }
   };
//
// drawing update
//
void face(float x0,float y0,float z0,float x1,float y1,float z1,
   float x2,float y2,float z2, float r, float g, float b, int *index) {
   v.data[(*index)++] = x0;
   v.data[(*index)++] = y0;
   v.data[(*index)++] = z0;
   v.data[(*index)++] = r;
   v.data[(*index)++] = g;
   v.data[(*index)++] = b;
   v.data[(*index)++] = x1;
   v.data[(*index)++] = y1;
   v.data[(*index)++] = z1;
   v.data[(*index)++] = r;
   v.data[(*index)++] = g;
   v.data[(*index)++] = b;
   v.data[(*index)++] = x2;
   v.data[(*index)++] = y2;
   v.data[(*index)++] = z2;
   v.data[(*index)++] = r;
   v.data[(*index)++] = g;
   v.data[(*index)++] = b;
   }
void update() {
   int index = 0;
   float scale = 2/(v.xmax-v.xmin);
   float d = v.size*v.d*scale;
   float xmid = (v.xmax+v.xmin)/2.0;
   float ymid = (v.ymax+v.ymin)/2.0;
   float zmid = (v.zmax+v.zmin)/2.0;
   float r,g,b;
   for (int i = 0; i < v.np; ++i) {
      if ((v.type[i] == 'b') || (v.type[i] == 't')) {
         r = 0.6;
         g = 0.6;
         b = 1.0;
         }
      else {
         r = 0.6+v.scale*v.s[i];
         g = 1.0-v.scale*v.s[i];
         b = 0.6;
         }
      float x = (v.px[i]-xmid)*scale;
      float y = (v.py[i]-ymid)*scale;
      float z = (zmid-v.pz[i])*scale; // ?
      face(x-d/2,y-d/2,z,x+d/2,y-d/2,z,x,y,z+d/2,r,g,b,&index);
      face(x-d/2,y+d/2,z,x,y,z+d/2,x+d/2,y+d/2,z,r,g,b,&index);
      face(x+d/2,y-d/2,z,x+d/2,y+d/2,z,x,y,z+d/2,r,g,b,&index);
      face(x-d/2,y+d/2,z,x-d/2,y-d/2,z,x,y,z+d/2,r,g,b,&index);
      face(x+d/2,y-d/2,z,x-d/2,y-d/2,z,x,y,z-d/2,r,g,b,&index);
      face(x-d/2,y+d/2,z,x+d/2,y+d/2,z,x,y,z-d/2,r,g,b,&index);
      face(x+d/2,y+d/2,z,x+d/2,y-d/2,z,x,y,z-d/2,r,g,b,&index);
      face(x-d/2,y-d/2,z,x-d/2,y+d/2,z,x,y,z-d/2,r,g,b,&index);
      }
   }
//
// main
//
int main(int argc, char** argv) {
   if (argc == 1) {
      cerr << "command line: strain source | strain_GL [arguments]" << endl;
      cerr << "   size=(particle size, 1=lattice pitch)" << endl;
      cerr << "   scale=(strain scale)" << endl;
      cerr << "   rx=(view x rotation radians)" << endl;
      cerr << "   ry=(view y rotation radians)" << endl;
      cerr << "   rz=(view z rotation radians)" << endl;
      cerr << "   sxyz=(view scale)" << endl;
      cerr << "   perspective=(3D projection, 0=orthographic)" << endl;
      return 1;
      }
   for (int i = 1; i < argc; ++i) {
      if (strstr(argv[i],"size=") != nullptr)
         v.size = stof(argv[i]+5);
      else if (strstr(argv[i],"scale=") != nullptr)
         v.scale = stof(argv[i]+6);
      else if (strstr(argv[i],"rx=") != nullptr)
         v.rx = stof(argv[i]+3);
      else if (strstr(argv[i],"ry=") != nullptr)
         v.ry = stof(argv[i]+3);
      else if (strstr(argv[i],"rz=") != nullptr)
         v.rz = stof(argv[i]+3);
      else if (strstr(argv[i],"sxyz=") != nullptr)
         v.sxyz = stof(argv[i]+5);
      else if (strstr(argv[i],"perspective=") != nullptr)
         v.perspective = stof(argv[i]+12);
      else {
         cerr << "strain_GL argument not recognized" << endl;
         return 1;
         }
      }
   //
   // read and parse inputs
   //
   string s;
   while (1) {
      getline(cin,s);
      if (s.compare("number:") == 0) {
         getline(cin,s);
         v.np = stoi(s);
         v.px = new float[v.np];
         v.py = new float[v.np];
         v.pz = new float[v.np];
         v.s = new float[v.np];
         v.type = new char[v.np];
         v.data = new float[CPP*v.np];
         }
      else if (s.compare("pitch:") == 0) {
         getline(cin,s);
         v.d = stof(s);
         }
      else if (s.compare("xmin:") == 0) {
         getline(cin,s);
         v.xmin = stof(s);
         }
      else if (s.compare("xmax:") == 0) {
         getline(cin,s);
         v.xmax = stof(s);
         }
      else if (s.compare("ymin:") == 0) {
         getline(cin,s);
         v.ymin = stof(s);
         }
      else if (s.compare("ymax:") == 0) {
         getline(cin,s);
         v.ymax = stof(s);
         }
      else if (s.compare("zmin:") == 0) {
         getline(cin,s);
         v.zmin = stof(s);
         }
      else if (s.compare("zmax:") == 0) {
         getline(cin,s);
         v.zmax = stof(s);
         }
      else if (s.compare("particles:") == 0) {
         for (int i = 0; i < v.np ; ++i) {
            cin.read((char *)&v.px[i],4);
            cin.read((char *)&v.py[i],4);
            cin.read((char *)&v.pz[i],4);
            }
         }
      else if (s.compare("type:") == 0) {
         for (int i = 0; i < v.np ; ++i)
            cin.read(&v.type[i],1);
         }
      else if (s.compare("continue:") == 0)
         break;
      }
   cerr << "strain_GL:" << endl;
   cerr << "   input " << v.np << " particles" << endl;
   //
   // initialize GL
   //
   GL g(v.width,v.height,(char *)"strain_GL");
   //
   // start main loop
   //
   int count = 0;
   char str[100];
   while (1) {
      getline(cin,s);
      if (s.compare("particles:") == 0)
         for (int i = 0; i < v.np ; ++i) {
            cin.read((char *)&v.px[i],4);
            cin.read((char *)&v.py[i],4);
            cin.read((char *)&v.pz[i],4);
            }
      else if (s.compare("strain:") == 0)
         for (int i = 0; i < v.np ; ++i)
            cin.read((char *)&v.s[i],4);
      else if (s.compare("continue:") == 0) {
         //
         // update GL
         //
         g.draw(update);
         //
         // uncomment to save images
         //
         //sprintf(str,"import -window strain_GL images/%04d.jpg",count++);
         //int result = system(str);
         }
      }
   }