diff options
Diffstat (limited to 'extern/mantaflow/helper/util/integrator.h')
-rw-r--r-- | extern/mantaflow/helper/util/integrator.h | 79 |
1 files changed, 79 insertions, 0 deletions
diff --git a/extern/mantaflow/helper/util/integrator.h b/extern/mantaflow/helper/util/integrator.h new file mode 100644 index 00000000000..5b1b02a5197 --- /dev/null +++ b/extern/mantaflow/helper/util/integrator.h @@ -0,0 +1,79 @@ +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2011 Tobias Pfaff, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * 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 Manta + +#endif |