integrator.h

Go to the documentation of this file.

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey 
 *
 * This program is free software, distributed under the terms of the
 * GNU General Public License (GPL) 
 * http://www.gnu.org/licenses
 *
 * Helper functions for simple integration
 *
 ******************************************************************************/

#ifndef _INTEGRATE_H
#define _INTEGRATE_H

#include <vector>
#include "vectorbase.h"
#include "kernel.h"

namespace Manta {
    
enum IntegrationMode { IntEuler=0, IntRK2, IntRK4 };
    
template<class VelKernel>
void integratePointSet(VelKernel& k, int mode) {
    typedef typename VelKernel::type0 PosType;
    PosType& x = k.getArg0();
    const std::vector<Vec3>& u = k.getRet();
    const int N = x.size();
        
    if (mode == IntEuler) {
        for(int i=0; i<N; i++) x[i].pos += u[i];
    } 
    else if (mode == IntRK2) {
        PosType x0(x);
        
        for(int i=0; i<N; i++) x[i].pos = x0[i].pos + 0.5*u[i];
        
        k.run();
        for(int i=0; i<N; i++) x[i].pos = x0[i].pos + u[i];
    } 
    else if (mode == IntRK4) {
        PosType x0(x);
        std::vector<Vec3> uTotal(u);
        
        for(int i=0; i<N; i++) x[i].pos = x0[i].pos + 0.5*u[i];
        
        k.run();
        for(int i=0; i<N; i++) {
            x[i].pos = x0[i].pos + 0.5*u[i];
            uTotal[i] += 2*u[i];
        }
        
        k.run();
        for(int i=0; i<N; i++) {
            x[i].pos = x0[i].pos + u[i];
            uTotal[i] += 2*u[i];
        }
        
        k.run();
        for(int i=0; i<N; i++) x[i].pos = x0[i].pos + (Real)(1./6.) * (uTotal[i] + u[i]);
    }
    else 
        errMsg("unknown integration type");
    
    //for(int i=0; i<N; i++) std::cout << x[i].pos.y-x[0].pos.y << std::endl;
    //std::cout << "<><><>" << std::endl;
}


} // namespace

#endif