Gradient Augmented Levelset Implementation in CPU & GPU
hermite-3d.cc (Latest change: Author:Lakshman Anumolu <acrlakshman@yahoo.co.in>, 2019-07-27 14:14:10 -0500, [commit: 0b0a8dc])
Go to the documentation of this file.
1 // Copyright 2019 Lakshman Anumolu, Raunak Bardia.
3 //
4 // Redistribution and use in source and binary forms, with or without
5 // modification, are permitted provided that the following conditions are
6 // met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the copyright holder nor the names of its contributors
16 // may be used to endorse or promote products derived from this software without
17 // specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
23 // HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
24 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
25 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
26 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
27 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
28 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
29 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 
32 #include "gals/cpu/interpolation/hermite.h"
33 
34 #include <math.h>
35 
36 template <typename T>
37 GALS::CPU::InterpolatedFields<GALS::CPU::Vec3<T>> GALS::INTERPOLATION::Hermite<T, GALS::CPU::Grid<T, 3>>::interpolate(
38  const GALS::CPU::Grid<typename GALS::CPU::Grid<T, 3>::value_type, GALS::CPU::Grid<T, 3>::dim>& grid,
39  const typename GALS::CPU::Grid<T, 3>::position_type& x_interp,
40  const GALS::CPU::Levelset<GALS::CPU::Grid<T, 3>, T>& levelset, const bool use_gradient_limiting)
41 {
42  GALS::CPU::InterpolatedFields<GALS::CPU::Vec3<T>> hermite_fields;
43 
44  using T_GRID = GALS::CPU::Grid<T, 3>;
45  const int dim = T_GRID::dim;
46  const int axis_x = 0, axis_y = 1, axis_z = 2;
47  const auto& axis_vectors = GALS::CPU::Grid<typename T_GRID::value_type, T_GRID::dim>::axis_vectors;
48  const GALS::CPU::Vec3<int> base_i_j_k = grid.baseNodeId(x_interp);
49  const GALS::CPU::Vec3<int> base_ip1_j_k =
50  GALS::CPU::Vec3<int>(base_i_j_k[0] + axis_vectors(0, 0), base_i_j_k[1], base_i_j_k[2]);
51  const GALS::CPU::Vec3<int> base_i_jp1_k =
52  GALS::CPU::Vec3<int>(base_i_j_k[0], base_i_j_k[1] + axis_vectors(1, 1), base_i_j_k[2]);
53  const GALS::CPU::Vec3<int> base_i_j_kp1 =
54  GALS::CPU::Vec3<int>(base_i_j_k[0], base_i_j_k[1], base_i_j_k[2] + axis_vectors(2, 2));
55  const GALS::CPU::Vec3<int> base_ip1_jp1_k =
56  GALS::CPU::Vec3<int>(base_i_j_k[0] + axis_vectors(0, 0), base_i_j_k[1] + axis_vectors(1, 1), base_i_j_k[2]);
57  const GALS::CPU::Vec3<int> base_i_jp1_kp1 =
58  GALS::CPU::Vec3<int>(base_i_j_k[0], base_i_j_k[1] + axis_vectors(1, 1), base_i_j_k[2] + axis_vectors(2, 2));
59  const GALS::CPU::Vec3<int> base_ip1_j_kp1 =
60  GALS::CPU::Vec3<int>(base_i_j_k[0] + axis_vectors(0, 0), base_i_j_k[1], base_i_j_k[2] + axis_vectors(2, 2));
61  const GALS::CPU::Vec3<int> base_ip1_jp1_kp1 = GALS::CPU::Vec3<int>(
62  base_i_j_k[0] + axis_vectors(0, 0), base_i_j_k[1] + axis_vectors(1, 1), base_i_j_k[2] + axis_vectors(2, 2));
63 
64  const typename T_GRID::position_type x_base = grid(base_i_j_k);
65  const auto& dx = grid.dX();
66  const auto& one_over_dx = grid.oneOverDX();
67 
68  GALS::CPU::Vec3<T> eta = (x_interp - x_base) * one_over_dx;
69 
70  const ControlPoints<T>& control_points_x = GALS::INTERPOLATION::get_control_points(
71  levelset.phiPrev()(base_i_j_k), levelset.psiPrev()(base_i_j_k)[axis_x], levelset.phiPrev()(base_ip1_j_k),
72  levelset.psiPrev()(base_ip1_j_k)[axis_x], dx[axis_x], use_gradient_limiting);
73  const ControlPoints<T>& control_points_y = GALS::INTERPOLATION::get_control_points(
74  levelset.phiPrev()(base_i_j_k), levelset.psiPrev()(base_i_j_k)[axis_y], levelset.phiPrev()(base_i_jp1_k),
75  levelset.psiPrev()(base_i_jp1_k)[axis_y], dx[axis_y], use_gradient_limiting);
76  const ControlPoints<T>& control_points_z = GALS::INTERPOLATION::get_control_points(
77  levelset.phiPrev()(base_i_j_k), levelset.psiPrev()(base_i_j_k)[axis_z], levelset.phiPrev()(base_i_j_kp1),
78  levelset.psiPrev()(base_i_j_kp1)[axis_z], dx[axis_z], use_gradient_limiting);
79  const ControlPoints<T>& control_points_xy = GALS::INTERPOLATION::get_control_points(
80  levelset.phiPrev()(base_ip1_j_k), levelset.psiPrev()(base_ip1_j_k)[axis_y], levelset.phiPrev()(base_ip1_jp1_k),
81  levelset.psiPrev()(base_ip1_jp1_k)[axis_y], dx[axis_y], use_gradient_limiting);
82  const ControlPoints<T>& control_points_yx = GALS::INTERPOLATION::get_control_points(
83  levelset.phiPrev()(base_i_jp1_k), levelset.psiPrev()(base_i_jp1_k)[axis_x], levelset.phiPrev()(base_ip1_jp1_k),
84  levelset.psiPrev()(base_ip1_jp1_k)[axis_x], dx[axis_x], use_gradient_limiting);
85  const ControlPoints<T>& control_points_xz = GALS::INTERPOLATION::get_control_points(
86  levelset.phiPrev()(base_ip1_j_k), levelset.psiPrev()(base_ip1_j_k)[axis_z], levelset.phiPrev()(base_ip1_j_kp1),
87  levelset.psiPrev()(base_ip1_j_kp1)[axis_z], dx[axis_z], use_gradient_limiting);
88  const ControlPoints<T>& control_points_zxy = GALS::INTERPOLATION::get_control_points(
89  levelset.phiPrev()(base_ip1_j_kp1), levelset.psiPrev()(base_ip1_j_kp1)[axis_y],
90  levelset.phiPrev()(base_ip1_jp1_kp1), levelset.psiPrev()(base_ip1_jp1_kp1)[axis_y], dx[axis_y],
91  use_gradient_limiting);
92  const ControlPoints<T>& control_points_xyz = GALS::INTERPOLATION::get_control_points(
93  levelset.phiPrev()(base_ip1_jp1_k), levelset.psiPrev()(base_ip1_jp1_k)[axis_z],
94  levelset.phiPrev()(base_ip1_jp1_kp1), levelset.psiPrev()(base_ip1_jp1_kp1)[axis_z], dx[axis_z],
95  use_gradient_limiting);
96  const ControlPoints<T>& control_points_zx = GALS::INTERPOLATION::get_control_points(
97  levelset.phiPrev()(base_i_j_kp1), levelset.psiPrev()(base_i_j_kp1)[axis_x], levelset.phiPrev()(base_ip1_j_kp1),
98  levelset.psiPrev()(base_ip1_j_kp1)[axis_x], dx[axis_x], use_gradient_limiting);
99  const ControlPoints<T>& control_points_zy = GALS::INTERPOLATION::get_control_points(
100  levelset.phiPrev()(base_i_j_kp1), levelset.psiPrev()(base_i_j_kp1)[axis_y], levelset.phiPrev()(base_i_jp1_kp1),
101  levelset.psiPrev()(base_i_jp1_kp1)[axis_y], dx[axis_y], use_gradient_limiting);
102  const ControlPoints<T>& control_points_yzx = GALS::INTERPOLATION::get_control_points(
103  levelset.phiPrev()(base_i_jp1_kp1), levelset.psiPrev()(base_i_jp1_kp1)[axis_x],
104  levelset.phiPrev()(base_ip1_jp1_kp1), levelset.psiPrev()(base_ip1_jp1_kp1)[axis_x], dx[axis_x],
105  use_gradient_limiting);
106  const ControlPoints<T>& control_points_yz = GALS::INTERPOLATION::get_control_points(
107  levelset.phiPrev()(base_i_jp1_k), levelset.psiPrev()(base_i_jp1_k)[axis_z], levelset.phiPrev()(base_i_jp1_kp1),
108  levelset.psiPrev()(base_i_jp1_kp1)[axis_z], dx[axis_z], use_gradient_limiting);
109  const ControlPoints<T>& control_points_y_z = GALS::INTERPOLATION::get_control_points(
110  levelset.psiPrev()(base_i_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_k)[1],
111  levelset.psiPrev()(base_i_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[1], dx[axis_y],
112  use_gradient_limiting);
113  const ControlPoints<T>& control_points_zy_z = GALS::INTERPOLATION::get_control_points(
114  levelset.psiPrev()(base_i_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[1],
115  levelset.psiPrev()(base_i_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[1], dx[axis_y],
116  use_gradient_limiting);
117  const ControlPoints<T>& control_points_z_x = GALS::INTERPOLATION::get_control_points(
118  levelset.psiPrev()(base_i_j_k)[axis_x], levelset.phiMixedDerivativesPrev()(base_i_j_k)[2],
119  levelset.psiPrev()(base_i_j_kp1)[axis_x], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[2], dx[axis_z],
120  use_gradient_limiting);
121  const ControlPoints<T>& control_points_x_y = GALS::INTERPOLATION::get_control_points(
122  levelset.psiPrev()(base_i_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_k)[0],
123  levelset.psiPrev()(base_ip1_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[0], dx[axis_x],
124  use_gradient_limiting);
125  const ControlPoints<T>& control_points_z_xy = GALS::INTERPOLATION::get_control_points(
126  levelset.phiMixedDerivativesPrev()(base_i_j_k)[0], levelset.phiMixedDerivativesPrev()(base_i_j_k)[3],
127  levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[0], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[3],
128  dx[axis_z], use_gradient_limiting);
129  const ControlPoints<T>& control_points_z_y = GALS::INTERPOLATION::get_control_points(
130  levelset.psiPrev()(base_i_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_k)[1],
131  levelset.psiPrev()(base_i_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[1], dx[axis_z],
132  use_gradient_limiting);
133  const ControlPoints<T>& control_points_x_z = GALS::INTERPOLATION::get_control_points(
134  levelset.psiPrev()(base_i_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_k)[2],
135  levelset.psiPrev()(base_ip1_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[2], dx[axis_x],
136  use_gradient_limiting);
137  const ControlPoints<T>& control_points_zx_z = GALS::INTERPOLATION::get_control_points(
138  levelset.psiPrev()(base_i_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[2],
139  levelset.psiPrev()(base_ip1_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[2], dx[axis_x],
140  use_gradient_limiting);
141  const ControlPoints<T>& control_points_zx_y = GALS::INTERPOLATION::get_control_points(
142  levelset.psiPrev()(base_i_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[0],
143  levelset.psiPrev()(base_ip1_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[0], dx[axis_x],
144  use_gradient_limiting);
145  const ControlPoints<T>& control_points_yx_y = GALS::INTERPOLATION::get_control_points(
146  levelset.psiPrev()(base_i_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[0],
147  levelset.psiPrev()(base_ip1_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[0], dx[axis_x],
148  use_gradient_limiting);
149  const ControlPoints<T>& control_points_yz_xy = GALS::INTERPOLATION::get_control_points(
150  levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[0], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[3],
151  levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[0], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[3],
152  dx[axis_z], use_gradient_limiting);
153  const ControlPoints<T>& control_points_yz_y = GALS::INTERPOLATION::get_control_points(
154  levelset.psiPrev()(base_i_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[1],
155  levelset.psiPrev()(base_i_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[1], dx[axis_z],
156  use_gradient_limiting);
157  const ControlPoints<T>& control_points_yx_z = GALS::INTERPOLATION::get_control_points(
158  levelset.psiPrev()(base_i_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[2],
159  levelset.psiPrev()(base_ip1_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[2], dx[axis_x],
160  use_gradient_limiting);
161  const ControlPoints<T>& control_points_yzx_y = GALS::INTERPOLATION::get_control_points(
162  levelset.psiPrev()(base_i_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[0],
163  levelset.psiPrev()(base_ip1_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[0], dx[axis_x],
164  use_gradient_limiting);
165  const ControlPoints<T>& control_points_yzx_z = GALS::INTERPOLATION::get_control_points(
166  levelset.psiPrev()(base_i_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[2],
167  levelset.psiPrev()(base_ip1_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[2], dx[axis_x],
168  use_gradient_limiting);
169  const ControlPoints<T>& control_points_xz_xy = GALS::INTERPOLATION::get_control_points(
170  levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[0], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[3],
171  levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[0], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[3],
172  dx[axis_z], use_gradient_limiting);
173  const ControlPoints<T>& control_points_xz_y = GALS::INTERPOLATION::get_control_points(
174  levelset.psiPrev()(base_ip1_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[1],
175  levelset.psiPrev()(base_ip1_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[1], dx[axis_z],
176  use_gradient_limiting);
177  const ControlPoints<T>& control_points_xyz_xy = GALS::INTERPOLATION::get_control_points(
178  levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[0], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[3],
179  levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[0], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[3],
180  dx[axis_z], use_gradient_limiting);
181  const ControlPoints<T>& control_points_xyz_y = GALS::INTERPOLATION::get_control_points(
182  levelset.psiPrev()(base_ip1_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[1],
183  levelset.psiPrev()(base_ip1_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[1], dx[axis_z],
184  use_gradient_limiting);
185  const ControlPoints<T>& control_points_xy_z = GALS::INTERPOLATION::get_control_points(
186  levelset.psiPrev()(base_ip1_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[1],
187  levelset.psiPrev()(base_ip1_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[1], dx[axis_y],
188  use_gradient_limiting);
189  const ControlPoints<T>& control_points_zxy_z = GALS::INTERPOLATION::get_control_points(
190  levelset.psiPrev()(base_ip1_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[1],
191  levelset.psiPrev()(base_ip1_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[1], dx[axis_y],
192  use_gradient_limiting);
193 
194  const T& c_21_x = control_points_x.c_21;
195  const T& c_12_x = control_points_x.c_12;
196  const T& c_21_y = control_points_y.c_21;
197  const T& c_12_y = control_points_y.c_12;
198  const T& c_30_z = control_points_z.c_30;
199  const T& c_21_z = control_points_z.c_21;
200  const T& c_12_z = control_points_z.c_12;
201  const T& c_03_z = control_points_z.c_03;
202  const T& c_21_xy = control_points_xy.c_21;
203  const T& c_12_xy = control_points_xy.c_12;
204  const T& c_21_yx = control_points_yx.c_21;
205  const T& c_12_yx = control_points_yx.c_12;
206  const T& c_30_xz = control_points_xz.c_30;
207  const T& c_21_xz = control_points_xz.c_21;
208  const T& c_12_xz = control_points_xz.c_12;
209  const T& c_03_xz = control_points_xz.c_03;
210  const T& c_21_zxy = control_points_zxy.c_21;
211  const T& c_12_zxy = control_points_zxy.c_12;
212  const T& c_30_xyz = control_points_xyz.c_30;
213  const T& c_21_xyz = control_points_xyz.c_21;
214  const T& c_12_xyz = control_points_xyz.c_12;
215  const T& c_03_xyz = control_points_xyz.c_03;
216  const T& c_21_zx = control_points_zx.c_21;
217  const T& c_12_zx = control_points_zx.c_12;
218  const T& c_21_zy = control_points_zy.c_21;
219  const T& c_12_zy = control_points_zy.c_12;
220  const T& c_21_yzx = control_points_yzx.c_21;
221  const T& c_12_yzx = control_points_yzx.c_12;
222  const T& c_30_yz = control_points_yz.c_30;
223  const T& c_21_yz = control_points_yz.c_21;
224  const T& c_12_yz = control_points_yz.c_12;
225  const T& c_03_yz = control_points_yz.c_03;
226  const T& c_21_y_z = control_points_y_z.c_21;
227  const T& c_12_y_z = control_points_y_z.c_12;
228  const T& c_21_zy_z = control_points_zy_z.c_21;
229  const T& c_12_zy_z = control_points_zy_z.c_12;
230  const T& c_21_z_x = control_points_z_x.c_21;
231  const T& c_12_z_x = control_points_z_x.c_12;
232  const T& c_21_x_y = control_points_x_y.c_21;
233  const T& c_12_x_y = control_points_x_y.c_12;
234  const T& c_21_z_xy = control_points_z_xy.c_21;
235  const T& c_12_z_xy = control_points_z_xy.c_12;
236  const T& c_21_z_y = control_points_z_y.c_21;
237  const T& c_12_z_y = control_points_z_y.c_12;
238  const T& c_21_x_z = control_points_x_z.c_21;
239  const T& c_12_x_z = control_points_x_z.c_12;
240  const T& c_21_zx_z = control_points_zx_z.c_21;
241  const T& c_12_zx_z = control_points_zx_z.c_12;
242  const T& c_21_zx_y = control_points_zx_y.c_21;
243  const T& c_12_zx_y = control_points_zx_y.c_12;
244  const T& c_21_yx_y = control_points_yx_y.c_21;
245  const T& c_12_yx_y = control_points_yx_y.c_12;
246  const T& c_21_yz_xy = control_points_yz_xy.c_21;
247  const T& c_12_yz_xy = control_points_yz_xy.c_12;
248  const T& c_21_yz_y = control_points_yz_y.c_21;
249  const T& c_12_yz_y = control_points_yz_y.c_12;
250  const T& c_21_yx_z = control_points_yx_z.c_21;
251  const T& c_12_yx_z = control_points_yx_z.c_12;
252  const T& c_21_yzx_y = control_points_yzx_y.c_21;
253  const T& c_12_yzx_y = control_points_yzx_y.c_12;
254  const T& c_21_yzx_z = control_points_yzx_z.c_21;
255  const T& c_12_yzx_z = control_points_yzx_z.c_12;
256  const T& c_21_xz_xy = control_points_xz_xy.c_21;
257  const T& c_12_xz_xy = control_points_xz_xy.c_12;
258  const T& c_21_xz_y = control_points_xz_y.c_21;
259  const T& c_12_xz_y = control_points_xz_y.c_12;
260  const T& c_21_xyz_xy = control_points_xyz_xy.c_21;
261  const T& c_12_xyz_xy = control_points_xyz_xy.c_12;
262  const T& c_21_xyz_y = control_points_xyz_y.c_21;
263  const T& c_12_xyz_y = control_points_xyz_y.c_12;
264  const T& c_21_xy_z = control_points_xy_z.c_21;
265  const T& c_12_xy_z = control_points_xy_z.c_12;
266  const T& c_21_zxy_z = control_points_zxy_z.c_21;
267  const T& c_12_zxy_z = control_points_zxy_z.c_12;
268 
269  // HERE
270  const T& hx_by_three = dx[axis_x] * (T)one_third;
271  const T& hy_by_three = dx[axis_y] * (T)one_third;
272  const T& hz_by_three = dx[axis_z] * (T)one_third;
273 
274  T bx[] = {B0(eta[0]), B1(eta[0]), B2(eta[0]), B3(eta[0])};
275  T by[] = {B0(eta[1]), B1(eta[1]), B2(eta[1]), B3(eta[1])};
276  T bz[] = {B0(eta[2]), B1(eta[2]), B2(eta[2]), B3(eta[2])};
277  T bx_prime[] = {B0_Prime(eta[0]), B1_Prime(eta[0]), B2_Prime(eta[0]), B3_Prime(eta[0])};
278  T by_prime[] = {B0_Prime(eta[1]), B1_Prime(eta[1]), B2_Prime(eta[1]), B3_Prime(eta[1])};
279  T bz_prime[] = {B0_Prime(eta[2]), B1_Prime(eta[2]), B2_Prime(eta[2]), B3_Prime(eta[2])};
280 
281  T control_points_all[] = {
282  c_30_z,
283  c_21_z,
284  c_12_z,
285  c_03_z,
286  c_21_y,
287  c_21_y + (dx[axis_z] * (T)one_third * c_21_y_z),
288  c_21_zy - (dx[axis_z] * (T)one_third * c_21_zy_z),
289  c_21_zy,
290  c_12_y,
291  c_12_y + (dx[axis_z] * (T)one_third * c_12_y_z),
292  c_12_zy - (dx[axis_z] * (T)one_third * c_12_zy_z),
293  c_12_zy,
294  c_30_yz,
295  c_21_yz,
296  c_12_yz,
297  c_03_yz,
298  c_21_x,
299  c_21_z + (dx[axis_x] * (T)one_third * c_21_z_x),
300  c_12_z + (dx[axis_x] * (T)one_third * c_12_z_x),
301  c_21_zx,
302  c_21_x + (dx[axis_y] * (T)one_third * c_21_x_y),
303  c_21_x + (hx_by_three * hy_by_three * c_21_z_xy) + (hy_by_three * c_21_z_y) + (hz_by_three * c_21_x_z),
304  c_21_zx + (hx_by_three * hy_by_three * c_12_z_xy) + (hy_by_three * c_12_z_y) - (hz_by_three * c_21_zx_z),
305  c_21_zx + (hy_by_three * c_21_zx_y),
306  c_21_yx - (hy_by_three * c_21_yx_y),
307  c_21_yx - (hx_by_three * hy_by_three * c_21_yz_xy) - (hy_by_three * c_21_yz_y) + (hz_by_three * c_21_yx_z),
308  c_21_yzx - (hx_by_three * hy_by_three * c_12_yz_xy) - (hy_by_three * c_12_yz_y) - (hz_by_three * c_21_yzx_z),
309  c_21_yzx - (hy_by_three * c_21_yzx_y),
310  c_21_yx,
311  c_21_yx + (hz_by_three * c_21_yx_z),
312  c_21_yzx - (hz_by_three * c_21_yzx_z),
313  c_21_yzx,
314  c_12_x,
315  c_12_x + (hz_by_three * c_12_x_z),
316  c_12_zx - (hz_by_three * c_12_zx_z),
317  c_12_zx,
318  c_12_x + (hy_by_three * c_12_x_y),
319  c_12_x - (hx_by_three * hy_by_three * c_21_xz_xy) + (hy_by_three * c_21_xz_y) + (hz_by_three * c_12_x_z),
320  c_12_zx - (hx_by_three * hy_by_three * c_12_xz_xy) + (hy_by_three * c_12_xz_y) - (hz_by_three * c_12_zx_z),
321  c_12_zx + (hy_by_three * c_12_zx_y),
322  c_12_yx - (hy_by_three * c_12_yx_y),
323  c_12_yx + (hx_by_three * hy_by_three * c_21_xyz_xy) - (hy_by_three * c_21_xyz_y) + (hz_by_three * c_12_yx_z),
324  c_12_yzx + (hx_by_three * hy_by_three * c_12_xyz_xy) - (hy_by_three * c_12_xyz_y) - (hz_by_three * c_12_yzx_z),
325  c_12_yzx - (hy_by_three * c_12_yzx_y),
326  c_12_yx,
327  c_12_yx + (hz_by_three * c_12_yx_z),
328  c_12_yzx - (hz_by_three * c_12_yzx_z),
329  c_12_yzx,
330  c_30_xz,
331  c_21_xz,
332  c_12_xz,
333  c_03_xz,
334  c_21_xy,
335  c_21_xy + (hz_by_three * c_21_xy_z),
336  c_21_zxy - (hz_by_three * c_21_zxy_z),
337  c_21_zxy,
338  c_12_xy,
339  c_12_xy + (hz_by_three * c_12_xy_z),
340  c_12_zxy - (hz_by_three * c_12_zxy_z),
341  c_12_zxy,
342  c_30_xyz,
343  c_21_xyz,
344  c_12_xyz,
345  c_03_xyz};
346 
347  for (int i = 0, l = 0; i < 4; ++i)
348  for (int j = 0; j < 4; ++j)
349  for (int k = 0; k < 4; ++k) {
350  hermite_fields.phi_interpolated += bx[i] * by[j] * bz[k] * control_points_all[l];
351  ++l;
352  }
353 
354  for (int i = 0, l = 0; i < 4; ++i)
355  for (int j = 0; j < 4; ++j)
356  for (int k = 0; k < 4; ++k) {
357  hermite_fields.psi_interpolated[axis_x] += bx_prime[i] * by[j] * bz[k] * control_points_all[l];
358  ++l;
359  }
360  hermite_fields.psi_interpolated[axis_x] *= one_over_dx[axis_x];
361 
362  for (int i = 0, l = 0; i < 4; ++i)
363  for (int j = 0; j < 4; ++j)
364  for (int k = 0; k < 4; ++k) {
365  hermite_fields.psi_interpolated[axis_y] += bx[i] * by_prime[j] * bz[k] * control_points_all[l];
366  ++l;
367  }
368  hermite_fields.psi_interpolated[axis_y] *= one_over_dx[axis_y];
369 
370  for (int i = 0, l = 0; i < 4; ++i)
371  for (int j = 0; j < 4; ++j)
372  for (int k = 0; k < 4; ++k) {
373  hermite_fields.psi_interpolated[axis_z] += bx[i] * by[j] * bz_prime[k] * control_points_all[l];
374  ++l;
375  }
376  hermite_fields.psi_interpolated[axis_z] *= one_over_dx[axis_z];
377 
378  return hermite_fields;
379 }
380 
381 template <typename T>
382 void GALS::INTERPOLATION::Hermite<T, GALS::CPU::Grid<T, 3>>::compute(
383  const GALS::CPU::Array<GALS::CPU::Grid<T, 3>, typename GALS::CPU::Grid<T, 3>::position_type>& x_interp,
384  GALS::CPU::Levelset<GALS::CPU::Grid<T, 3>, T>& levelset)
385 {
386  typedef GALS::CPU::Grid<T, 3> T_GRID;
387 
388  const GALS::CPU::Vec3<int> num_cells_interp = x_interp.numCells();
389  const T_GRID& grid = levelset.grid();
390  const GALS::CPU::Vec3<typename T_GRID::value_type> dx = grid.dX();
391  const auto& axis_vectors = GALS::CPU::Grid<typename T_GRID::value_type, T_GRID::dim>::axis_vectors;
392 
393  for (int i = 0; i < num_cells_interp[0]; ++i)
394  for (int j = 0; j < num_cells_interp[1]; ++j)
395  for (int k = 0; k < num_cells_interp[2]; ++k) {
396  const auto& hermite_fields = this->interpolate(grid, x_interp(i, j, k), levelset);
397 
398  levelset.phiInterpPrev()(i, j, k) = hermite_fields.phi_interpolated;
399  levelset.psiInterpPrev()(i, j, k) = hermite_fields.psi_interpolated;
400  }
401 }
402 
403 template class GALS::INTERPOLATION::Hermite<double, GALS::CPU::Grid<double, 3>>;
GALS::INTERPOLATION::B0
T B0(T eta)
cubic Bernstein polynomials
Definition: hermite.h:57
GALS::INTERPOLATION::B0_Prime
T B0_Prime(T eta)
Definition: hermite.h:81
GALS::CPU::Grid::dX
const Vec3< T > & dX() const
Definition: grid.cc:144
GALS::CPU::Grid::value_type
T value_type
Definition: grid.h:71
GALS::INTERPOLATION::B2_Prime
T B2_Prime(T eta)
Definition: hermite.h:93
GALS::one_third
static const double one_third
Definition: utilities.h:58
GALS::INTERPOLATION::get_control_points
static struct ControlPoints< T > get_control_points(const T &phi_0, const T &phi_x_0, const T &phi_1, const T &phi_x_1, const T &dx, bool use_gradient_limiting=false)
Definition: hermite.h:116
GALS::INTERPOLATION::B3_Prime
T B3_Prime(T eta)
Definition: hermite.h:99
GALS::INTERPOLATION::B1_Prime
T B1_Prime(T eta)
Definition: hermite.h:87