Changeset View
Changeset View
Standalone View
Standalone View
intern/mantaflow/intern/manta_pp/omp/util/integrator.h
- This file was added.
| /****************************************************************************** | |||||
| * | |||||
| * 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 }; | |||||
| //! Integrate a particle set with a given velocity kernel | |||||
| 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 | |||||