blob: 5b1b02a519757d88047e43ee91fb4b219fc261ba (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
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
|
