vlmain.cpp

#include <boost/program_options.hpp>
#include <boost/thread.hpp>
#include <boost/signal.hpp>
#include <boost/bind.hpp>
#include <boost/format.hpp>
namespace po = boost::program_options;

#include <vlCore/VisualizationLibrary.hpp>
#include <vlGLUT/GLUTWindow.hpp>
#include <vlCore/VLSettings.hpp>

#include "logger/logger.h"
#include "logger/stdiowriter.h"
#include "logger/syslogwriter.h"

#include "core/common.h"
#include "core/fields.hpp"
#include "core/scalars.hpp"
#include "core/macsmoke.hpp"
#include "core/inssmoke.hpp"
#include "core/pointrender.hpp"
#include "core/volumerender.hpp"
#include "render/glutapp.h"
#include "render/drawable.h"
#include "io/xmlparser.h"

#ifndef TARGET_VERSION_MAJOR
#define TARGET_VERSION_MAJOR 99
#endif

#ifndef TARGET_VERSION_MINOR
#define TARGET_VERSION_MINOR 99
#endif

static fdl::MACSmoke<float>* grid = NULL;
static fdl::INSSmoke<float>* solver = NULL;
fdl::Constant<float>* gravityField = NULL;
// static fdl::FieldBase<float>* advectField = NULL;
static fdl::ScalarBase<float>* smokeSrc = NULL;
static fdl::ScalarBase<float>* tempSrc = NULL;
static vl::ref<fdl::VolumeRenderer> applet;

template<class T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
{
        std::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " ")); 
        return os;
}

void simulateSmoke()
{
        for (int t = 0; t < 200; ++t) {
                
                INFO() << boost::format("T= %g, S= %g, p= %g, u= %g, v= %g, w= %g") \
                % grid->supP() % grid->supT() % grid->supS() % grid->supU() % grid->supV() % grid->supW();
                
                INFO() << "Advection step";
                solver->advect();
                
                INFO() << "Adding external force";
                solver->addfield();
                
                INFO() << "Adding source terms";
                solver->addsources();
                
                INFO() << "Pressure update step";
                solver->project();
        }
}

int main(int argc, char **argv)
{
        fdl::Logger out;
        // figure out how to automate this output...
        std::cout << "FDL version " << TARGET_VERSION_MAJOR << "." << TARGET_VERSION_MINOR << " of " << __DATE__ << " at " << __TIME__ << std::endl;
        
        fdl::SyslogWriter* syslogger = new fdl::SyslogWriter();
        fdl::StdOutWriter* stdlogger = new fdl::StdOutWriter();
        stdlogger->setFormat("%H:%M:%S");
        fdl::Logger::setIdentity("FDL");
        // fdl::Logger::RegisterWriter(syslogger);
        fdl::Logger::registerWriter(stdlogger);
        fdl::Logger::setLevel(fdl::Logger::WARN | fdl::Logger::ERROR);
        
        std::string output_prefix = "density_export_";
        
        int window_width, window_height;
        int use_glut_window = 1;
        std::string render_exporter = "";
        std::string levelLabel = "";
        std::string inputfile = "";
        int level = 0;
        
        float dx = 0.1;
        float dt = 0.1;
        float rho = 1.0;
        std::vector<int> grid_dims(3, 50);
        float solidTemp = 0.0;
        float smokeSrcTemp = 2.0;
        
        float gvec[3] = {0.0, 0.0, -1.0}; // gravitational field (constant)

        float cg_tol = fdl::EPSILON; 
        unsigned cg_maxiter = 200;
        float ode_tol = 1e-2;
        
        try {
                po::options_description desc("Allowed options");
                desc.add_options()
                ("help,H", "produce help message")
                ("input-file,C", po::value<std::string>(&inputfile), "input configuration file")
                ("output-format,O", po::value<std::string>(), "output format")
                ("output-name,N", po::value<std::string>(&output_prefix), "output file name PREFIX")
                ("win-width,W", po::value<int>(&window_width)->default_value(512), "Glut window width")
                ("win-height,E", po::value<int>(&window_height)->default_value(512), "Glut window height")
                ("use-glut,L", po::value<int>(&use_glut_window)->default_value(1), "Use Glut rendering window")
                ("verbose,V", po::value<int>(&level), "logging level");
                
                po::variables_map vm;
                po::store(po::command_line_parser(argc, argv).options(desc).style(po::command_line_style::default_style 
                        | po::command_line_style::allow_slash_for_short 
                        | po::command_line_style::allow_long_disguise).run(), vm);
                
                po::notify(vm);
                
                if (vm.count("help")) {
                        std::cout << "Usage: options_description [options]" << std::endl;
                        std::cout << desc << std::endl;
                        return 0;
                }
                
                if (vm.count("output-name")) {
                        std::cout << "Output name is: " << vm["output-name"].as<std::string>() << std::endl;
                }
                
                if (vm.count("output-format")) {
                        std::cout << "Output format is: " << vm["output-format"].as<std::string>() << std::endl;
                }
                
                if(level >= 0){
                        levelLabel = fdl::Logger::loggerLevelAsString(fdl::Logger::INFO);
                        fdl::Logger::pushLevel(fdl::Logger::INFO);
                }
                if(level >= 1){
                        levelLabel = fdl::Logger::loggerLevelAsString(fdl::Logger::DEV);
                        fdl::Logger::pushLevel(fdl::Logger::DEV);
                }
                if(level >= 2){
                        levelLabel = fdl::Logger::loggerLevelAsString(fdl::Logger::DEBUG);
                        fdl::Logger::pushLevel(fdl::Logger::DEBUG);
                }
                                
                if (vm.count("renderer")) {
                        std::cout << "Output name is: " << vm["output-name"].as<std::string>() << std::endl;
                }

                gravityField = new fdl::Constant<float>(3, gvec);

                float x_rel = 0.5, y_rel = 0.5;
                float sph_R0[3] = {grid_dims[0] * dx * x_rel, grid_dims[1] * dx * y_rel, dx};
                float sph_r = dx * 3;

                float cyl_XY[2] = {grid_dims[0] * dx * x_rel, grid_dims[1] * dx * y_rel};
                float cyl_z[2] = {dx, 3 * dx};
                float cyl_r = dx * 4;

                smokeSrc = new fdl::Cylindrical<float>(3, cyl_XY, cyl_r, cyl_z[0], cyl_z[1], 0.1 );
                tempSrc = new fdl::Gaussian<float>(3, sph_R0, sph_r, smokeSrcTemp);

                grid = new fdl::MACSmoke<float>(grid_dims[0], grid_dims[1], grid_dims[2], dx, solidTemp);
                // solver = new fdl::INSSmoke<float>(grid, gravityField, rho, dt, tempScalar, smokeSrc);
                
                
                INFO() << "Initializing for fluid at rest in solid container ";
                
                float X_tmp[3];
                for (int x = 1; x <= grid->xdim(); ++x) {
                        for (int y = 1; y <= grid->ydim(); ++y) {
                                for (int z = 1; z <= grid->zdim(); ++z) 
                                {
                                        // at the boundary...
                                        if (x == 1 || x == grid->xdim() ||
                                                 y == 1 || y == grid->ydim() ||
                                                 z == 1 || z == grid->zdim())
                                        {
                                                grid->operator()(x, y, z) = fdl::SOLID;
                                                grid->p(x, y, z) = 0.0;
                                                grid->u(x, y, z) = 0.0;
                                                grid->u(x+1, y, z) = 0.0;
                                                grid->v(x, y, z) = 0.0;
                                                grid->v(x, y+1, z) = 0.0;
                                                grid->w(x, y, z) = 0.0;
                                                grid->w(x, y, z+1) = 0.0;
                                                grid->T(x, y, z) = 1.0; // this is a hack
                                                grid->s(x, y, z) = 0.0;
                                        } else {
                                                grid->operator()(x, y, z) = fdl::FLUID;
                                                grid->p(x, y, z) = 1.0;
                                                grid->u(x, y, z) = 0.0;
                                                grid->v(x, y, z) = 0.0;
                                                grid->w(x, y, z) = 0.0;
                                                grid->coordsP(x, y, z, X_tmp);
                                                grid->T(x, y, z) = (*tempSrc)(0.0, X_tmp);
                                                // grid->T(x, y, z) = 0.0;
                                                grid->s(x, y, z) = (*smokeSrc)(0.0, X_tmp);
                                        }
                                }
                        }
                }
                
                solver = new fdl::INSSmoke<float>(grid, gravityField, rho, dt, tempSrc, smokeSrc);
                solver->setStepTolerance(ode_tol);
                solver->setStepRule(fdl::RK2);
                
                INFO() << "Logging level is " << levelLabel << " (" << level << ")";
                INFO() << "Grid dimensions are: " << grid->xdim() << " x " << grid->ydim() << " x " << grid->zdim();
                INFO() << "Cell width is: " << grid->dx();
                INFO() << "Time step is:  " << solver->dt();
                INFO() << "CG solver tolerance is: " << solver->getCGTolerance();
                INFO() << "CG solver max iterations is: " << solver->getCGMaxIter();
                INFO() << "Fluid density coeff. is:  " << solver->rho();
                
                if(use_glut_window>0){
                        boost::thread t(simulateSmoke);
                        
                        int pargc = argc;
                        glutInit( &pargc, argv );
                
                        // init Visualization Library
                        vl::VisualizationLibrary::init();
                        vl::globalSettings()->setDataPath(vl::String("."));
                
                        // setup the OpenGL context format
                        vl::OpenGLContextFormat format;
                        format.setDoubleBuffer(true);
                        format.setRGBABits( 8,8,8,8 );
                        format.setDepthBufferBits(24);
                        format.setStencilBufferBits(8);
                        format.setFullscreen(false);
                        format.setMultisampleSamples(8);
                        format.setMultisample(true);
                                        
                        // create the applet to be run
                        //applet = new fdl::PointRenderer();
                        applet = new fdl::VolumeRenderer();
                        applet->initialize();
                        applet->setGrid(grid);
                        //applet->setSolver(solver);
                        // create a native GLUT window
                        vl::ref<vlGLUT::GLUTWindow> glut_window = new vlGLUT::GLUTWindow;
                        // bind the applet so it receives all the GUI events related to the OpenGLContext
                        glut_window->addEventListener(applet.get());            
                        // target the window so we can render on it
                        applet->rendering()->as<Rendering>()->renderer()->setRenderTarget( glut_window->renderTarget() );
                        // black background
                        applet->rendering()->as<Rendering>()->camera()->viewport()->setClearColor( vl::black );
                        // define the camera position and orientation
                        vl::vec3 eye    = vl::vec3(0,30,0);     // camera position
                        vl::vec3 target = vl::vec3(0,0,0);      // point the camera is looking at
                        vl::vec3 up     = vl::vec3(0,1,0);      // up direction
                        
                        // --- trackball manipulator ---
                        vl::ref<vl::TrackballManipulator> trackball = new vl::TrackballManipulator;
                        // bind the camera to be manipulated
                        trackball->setCamera( applet->rendering()->as<Rendering>()->camera() );
                        // manipulate the camera's transform
                        trackball->setTransform(NULL);
                        // rotation pivot coordinates
                        trackball->setPivot( vl::vec3(0,0,0) );
                        // install the trackball
                        glut_window->addEventListener(trackball.get());
                        
                        AABB volume_box( vec3( 0,0,0 ), vec3( 1,1,1 ) );
                        trackball->adjustView(volume_box, vl::vec3(10,10,10), up);
                
                        // Initialize the OpenGL context and window properties
                        glut_window->initGLUTWindow( "FDL", format, 0, 0, window_width, window_height );
                
                        glutMainLoop();
                }
                else{
                        simulateSmoke();
                }
                
        }
        catch(std::exception& e) {
                ERROR() << " * error: " << e.what();
                return 1;
        }
        catch(...) {
                ERROR() << "Exception of unknown type!";
        }
        
        /* 
         Idea: move this stuff into the class constructors 
        
        std::cout << "Grid dimensions are: " << grid_dims[0] << " x " << grid_dims[1] << " x " << grid_dims[2] << std::endl;
        std::cout << "Cell width is: " << dx << std::endl;
        // and so on...
         */
        
        delete solver;
        delete grid;
                
        return 0;
}