-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathnumerics.hpp
143 lines (120 loc) · 5.71 KB
/
numerics.hpp
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
/*
* numerics.hpp
*
* Numerical building blocks that operate on Field() objects.
* Also defines useful typedefs to avoid template usage. IE DenseScalarField2 is a 2D non-sparse scalar field.
*
* Example:
* DenseVectorField3 grad_of_x;
* grad_of_x.createEntireLevel(8);
* grad<3>(x); // x is instance of ScalarField3
*
* Created on: Dec 27, 2016
* Author: Christian Huettig
*/
#pragma once
typedef HCS<1> H1; // 1D
typedef HCS<2> H2; // 2D
typedef HCS<3> H3; // 3D
typedef HCS<4> H4; // 4D
typedef HCS<5> H5; // 5D
typedef Tensor1<data_t, 2> Vec2; // Single "vector" in 2D / 3D. Similar to old Point<T>
typedef Tensor1<data_t, 3> Vec3;
typedef Tensor1<data_t, 4> Vec4;
typedef Tensor1<data_t, 5> Vec5;
typedef Field<data_t, H1> ScalarField1; // 1D scalar field type
typedef Field<data_t, H2> ScalarField2;
typedef Field<data_t, H3> ScalarField3;
typedef Field<data_t, H4> ScalarField4;
typedef Field<data_t, H5> ScalarField5;
template <unsigned dimensions> using ScalarField = Field<data_t, HCS<dimensions> >;
typedef Field<data_t, H1> VectorField1; // 1D scalar field type
typedef Field<Vec2, H2> VectorField2;
typedef Field<Vec3, H3> VectorField3;
typedef Field<Vec4, H4> VectorField4;
typedef Field<Vec5, H5> VectorField5;
template <unsigned dimensions> using VectorField = Field<Tensor1<data_t, dimensions>, HCS<dimensions> >;
typedef DenseField<data_t, H1> DenseScalarField1; // 1D scalar field type
typedef DenseField<data_t, H2> DenseScalarField2;
typedef DenseField<data_t, H3> DenseScalarField3;
typedef DenseField<data_t, H4> DenseScalarField4;
typedef DenseField<data_t, H5> DenseScalarField5;
template <unsigned dimensions> using DenseScalarField = DenseField<data_t, HCS<dimensions> >;
typedef DenseField<data_t, H1> DenseVectorField1; // 1D Vector field type
typedef DenseField<Vec2, H2> DenseVectorField2;
typedef DenseField<Vec3, H3> DenseVectorField3;
typedef DenseField<Vec4, H4> DenseVectorField4;
typedef DenseField<Vec5, H5> DenseVectorField5;
template <unsigned dimensions> using DenseVectorField = DenseField<Tensor1<data_t, dimensions>, HCS<dimensions> >;
typedef SparseField<data_t, H1> SparseScalarField1; // 1D scalar field type
typedef SparseField<data_t, H2> SparseScalarField2;
typedef SparseField<data_t, H3> SparseScalarField3;
typedef SparseField<data_t, H4> SparseScalarField4;
typedef SparseField<data_t, H5> SparseScalarField5;
template <unsigned dimensions> using SparseScalarField = SparseField<data_t, HCS<dimensions> >;
typedef SparseField<data_t, H1> SparseVectorField1; // 1D Vector field type
typedef SparseField<Vec2, H2> SparseVectorField2;
typedef SparseField<Vec3, H3> SparseVectorField3;
typedef SparseField<Vec4, H4> SparseVectorField4;
typedef SparseField<Vec5, H5> SparseVectorField5;
template <unsigned dimensions> using SparseVectorField = SparseField<Tensor1<data_t, dimensions>, HCS<dimensions> >;
// Turn 2 HCS coords into a vector containing pos(a)-pos(b)
// HCS must be provided as it may contain non-default scaling info
template <unsigned dimensions>
Tensor1<data_t, dimensions> c2v(HCS<dimensions> &hcs, coord_t a, coord_t b) {
Tensor1<data_t, dimensions> pos_a(hcs.getPosition(a));
Tensor1<data_t, dimensions> pos_b(hcs.getPosition(b));
return pos_a - pos_b;
}
// Dimension-independent gradient operator
// Takes a scalar field as input and returns the vector field
// Example: VectorField3 grad_of_x = grad<3>(x); // x is instance of ScalarField3
template<int dimension>
void grad(Field<data_t, HCS<dimension> > &source, Field<Tensor1<data_t, dimension>, HCS<dimension> > &result) {
typedef Tensor1<data_t, dimension> Vec;
typedef HCS<dimension> HX;
typedef Field<Vec, HX> VectorField;
typedef Field<data_t, HX> ScalarField;
result.createEntireLevel(source.getHighestLevel());
// strange calling convention because template depends now on another template
result.template convert<data_t>(source, [](coord_t c, ScalarField &source)->Vec {
HX &h = source.hcs;
Vec gradient;
// Gradient calculation with finite-difference stencil
level_t l = h.GetLevel(c);
data_t dist = 4 * (h.scales[0] / data_t(1U << l)); // neighbor distance at that level, assuming all scales equal, * 2 because gradient is 2nd order
for (int n_idx = 0; n_idx < h.parts; n_idx++) { // Traverse all neighbors
coord_t c_ne = h.getNeighbor(c, n_idx);
data_t ne_val = source.get(c_ne); // will respect boundary condition
gradient[n_idx >> 1] += (n_idx & 1 ? -ne_val : ne_val) / dist;
}
return gradient;
});
//result.propagate();
}
// Dimension-independent divergence operator
// Takes a scalar field as input and returns the vector field
// Example: VectorField3 grad_of_x = grad<3>(x); // x is instance of ScalarField3
template<int dimension>
void div(Field<Tensor1<data_t, dimension>, HCS<dimension> > &source, Field<data_t, HCS<dimension> > &result) {
typedef Tensor1<data_t, dimension> Vec;
typedef HCS<dimension> HX;
typedef Field<data_t, HX> ScalarField;
typedef Field<Vec, HX> VectorField;
result.createEntireLevel(source.getHighestLevel());
// strange calling convention because template depends now on another template
result.template convert<Vec>(source, [](coord_t c, VectorField &source)->data_t {
HX &h = source.hcs;
data_t divergence = 0;
// Gradient calculation with finite-difference stencil
level_t l = h.GetLevel(c);
data_t dist = 4 * (h.scales[0] / data_t(1U << l)); // neighbor distance at that level, assuming all scales equal, * 2 because gradient is 2nd order
for (int n_idx = 0; n_idx < h.parts; n_idx++) { // Traverse all neighbors
coord_t c_ne = h.getNeighbor(c, n_idx);
data_t ne_val = source.get(c_ne)[n_idx >> 1]; // will respect boundary condition
divergence += (n_idx & 1 ? -ne_val : ne_val) / dist;
}
return divergence;
});
//result.propagate();
}