Line data Source code
1 : ///////////////////////////////////////////////////////////////////////////////
2 : // 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.
30 : ///////////////////////////////////////////////////////////////////////////////
31 :
32 : #include "gals/cpu/interpolation/hermite.h"
33 :
34 : #include <math.h>
35 :
36 : template <typename T>
37 4993 : 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 4993 : 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 9986 : const GALS::CPU::Vec3<int> base_i_j_k = grid.baseNodeId(x_interp);
49 14979 : const GALS::CPU::Vec3<int> base_ip1_j_k =
50 4993 : GALS::CPU::Vec3<int>(base_i_j_k[0] + axis_vectors(0, 0), base_i_j_k[1], base_i_j_k[2]);
51 14979 : const GALS::CPU::Vec3<int> base_i_jp1_k =
52 4993 : GALS::CPU::Vec3<int>(base_i_j_k[0], base_i_j_k[1] + axis_vectors(1, 1), base_i_j_k[2]);
53 9986 : const GALS::CPU::Vec3<int> base_i_j_kp1 =
54 4993 : GALS::CPU::Vec3<int>(base_i_j_k[0], base_i_j_k[1], base_i_j_k[2] + axis_vectors(2, 2));
55 19972 : const GALS::CPU::Vec3<int> base_ip1_jp1_k =
56 9986 : 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 14979 : const GALS::CPU::Vec3<int> base_i_jp1_kp1 =
58 9986 : 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 14979 : const GALS::CPU::Vec3<int> base_ip1_j_kp1 =
60 9986 : 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 19972 : const GALS::CPU::Vec3<int> base_ip1_jp1_kp1 = GALS::CPU::Vec3<int>(
62 14979 : 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 9986 : const typename T_GRID::position_type x_base = grid(base_i_j_k);
65 4993 : const auto& dx = grid.dX();
66 4993 : const auto& one_over_dx = grid.oneOverDX();
67 :
68 9986 : GALS::CPU::Vec3<T> eta = (x_interp - x_base) * one_over_dx;
69 :
70 19972 : const ControlPoints<T>& control_points_x = GALS::INTERPOLATION::get_control_points(
71 9986 : levelset.phiPrev()(base_i_j_k), levelset.psiPrev()(base_i_j_k)[axis_x], levelset.phiPrev()(base_ip1_j_k),
72 14979 : levelset.psiPrev()(base_ip1_j_k)[axis_x], dx[axis_x], use_gradient_limiting);
73 14979 : const ControlPoints<T>& control_points_y = GALS::INTERPOLATION::get_control_points(
74 9986 : levelset.phiPrev()(base_i_j_k), levelset.psiPrev()(base_i_j_k)[axis_y], levelset.phiPrev()(base_i_jp1_k),
75 9986 : levelset.psiPrev()(base_i_jp1_k)[axis_y], dx[axis_y], use_gradient_limiting);
76 14979 : const ControlPoints<T>& control_points_z = GALS::INTERPOLATION::get_control_points(
77 9986 : levelset.phiPrev()(base_i_j_k), levelset.psiPrev()(base_i_j_k)[axis_z], levelset.phiPrev()(base_i_j_kp1),
78 9986 : levelset.psiPrev()(base_i_j_kp1)[axis_z], dx[axis_z], use_gradient_limiting);
79 14979 : const ControlPoints<T>& control_points_xy = GALS::INTERPOLATION::get_control_points(
80 9986 : levelset.phiPrev()(base_ip1_j_k), levelset.psiPrev()(base_ip1_j_k)[axis_y], levelset.phiPrev()(base_ip1_jp1_k),
81 9986 : levelset.psiPrev()(base_ip1_jp1_k)[axis_y], dx[axis_y], use_gradient_limiting);
82 14979 : const ControlPoints<T>& control_points_yx = GALS::INTERPOLATION::get_control_points(
83 9986 : levelset.phiPrev()(base_i_jp1_k), levelset.psiPrev()(base_i_jp1_k)[axis_x], levelset.phiPrev()(base_ip1_jp1_k),
84 9986 : levelset.psiPrev()(base_ip1_jp1_k)[axis_x], dx[axis_x], use_gradient_limiting);
85 14979 : const ControlPoints<T>& control_points_xz = GALS::INTERPOLATION::get_control_points(
86 9986 : levelset.phiPrev()(base_ip1_j_k), levelset.psiPrev()(base_ip1_j_k)[axis_z], levelset.phiPrev()(base_ip1_j_kp1),
87 9986 : levelset.psiPrev()(base_ip1_j_kp1)[axis_z], dx[axis_z], use_gradient_limiting);
88 14979 : const ControlPoints<T>& control_points_zxy = GALS::INTERPOLATION::get_control_points(
89 9986 : levelset.phiPrev()(base_ip1_j_kp1), levelset.psiPrev()(base_ip1_j_kp1)[axis_y],
90 9986 : levelset.phiPrev()(base_ip1_jp1_kp1), levelset.psiPrev()(base_ip1_jp1_kp1)[axis_y], dx[axis_y],
91 : use_gradient_limiting);
92 14979 : const ControlPoints<T>& control_points_xyz = GALS::INTERPOLATION::get_control_points(
93 9986 : levelset.phiPrev()(base_ip1_jp1_k), levelset.psiPrev()(base_ip1_jp1_k)[axis_z],
94 9986 : levelset.phiPrev()(base_ip1_jp1_kp1), levelset.psiPrev()(base_ip1_jp1_kp1)[axis_z], dx[axis_z],
95 : use_gradient_limiting);
96 14979 : const ControlPoints<T>& control_points_zx = GALS::INTERPOLATION::get_control_points(
97 9986 : levelset.phiPrev()(base_i_j_kp1), levelset.psiPrev()(base_i_j_kp1)[axis_x], levelset.phiPrev()(base_ip1_j_kp1),
98 9986 : levelset.psiPrev()(base_ip1_j_kp1)[axis_x], dx[axis_x], use_gradient_limiting);
99 14979 : const ControlPoints<T>& control_points_zy = GALS::INTERPOLATION::get_control_points(
100 9986 : levelset.phiPrev()(base_i_j_kp1), levelset.psiPrev()(base_i_j_kp1)[axis_y], levelset.phiPrev()(base_i_jp1_kp1),
101 9986 : levelset.psiPrev()(base_i_jp1_kp1)[axis_y], dx[axis_y], use_gradient_limiting);
102 14979 : const ControlPoints<T>& control_points_yzx = GALS::INTERPOLATION::get_control_points(
103 9986 : levelset.phiPrev()(base_i_jp1_kp1), levelset.psiPrev()(base_i_jp1_kp1)[axis_x],
104 9986 : levelset.phiPrev()(base_ip1_jp1_kp1), levelset.psiPrev()(base_ip1_jp1_kp1)[axis_x], dx[axis_x],
105 : use_gradient_limiting);
106 14979 : const ControlPoints<T>& control_points_yz = GALS::INTERPOLATION::get_control_points(
107 9986 : levelset.phiPrev()(base_i_jp1_k), levelset.psiPrev()(base_i_jp1_k)[axis_z], levelset.phiPrev()(base_i_jp1_kp1),
108 9986 : levelset.psiPrev()(base_i_jp1_kp1)[axis_z], dx[axis_z], use_gradient_limiting);
109 4993 : const ControlPoints<T>& control_points_y_z = GALS::INTERPOLATION::get_control_points(
110 9986 : levelset.psiPrev()(base_i_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_k)[1],
111 14979 : levelset.psiPrev()(base_i_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[1], dx[axis_y],
112 : use_gradient_limiting);
113 4993 : const ControlPoints<T>& control_points_zy_z = GALS::INTERPOLATION::get_control_points(
114 9986 : levelset.psiPrev()(base_i_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[1],
115 9986 : levelset.psiPrev()(base_i_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[1], dx[axis_y],
116 : use_gradient_limiting);
117 4993 : const ControlPoints<T>& control_points_z_x = GALS::INTERPOLATION::get_control_points(
118 9986 : levelset.psiPrev()(base_i_j_k)[axis_x], levelset.phiMixedDerivativesPrev()(base_i_j_k)[2],
119 9986 : levelset.psiPrev()(base_i_j_kp1)[axis_x], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[2], dx[axis_z],
120 : use_gradient_limiting);
121 4993 : const ControlPoints<T>& control_points_x_y = GALS::INTERPOLATION::get_control_points(
122 9986 : levelset.psiPrev()(base_i_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_k)[0],
123 9986 : levelset.psiPrev()(base_ip1_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[0], dx[axis_x],
124 : use_gradient_limiting);
125 4993 : const ControlPoints<T>& control_points_z_xy = GALS::INTERPOLATION::get_control_points(
126 9986 : levelset.phiMixedDerivativesPrev()(base_i_j_k)[0], levelset.phiMixedDerivativesPrev()(base_i_j_k)[3],
127 9986 : levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[0], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[3],
128 9986 : dx[axis_z], use_gradient_limiting);
129 4993 : const ControlPoints<T>& control_points_z_y = GALS::INTERPOLATION::get_control_points(
130 9986 : levelset.psiPrev()(base_i_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_k)[1],
131 9986 : levelset.psiPrev()(base_i_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[1], dx[axis_z],
132 : use_gradient_limiting);
133 4993 : const ControlPoints<T>& control_points_x_z = GALS::INTERPOLATION::get_control_points(
134 9986 : levelset.psiPrev()(base_i_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_k)[2],
135 9986 : levelset.psiPrev()(base_ip1_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[2], dx[axis_x],
136 : use_gradient_limiting);
137 4993 : const ControlPoints<T>& control_points_zx_z = GALS::INTERPOLATION::get_control_points(
138 9986 : levelset.psiPrev()(base_i_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[2],
139 9986 : levelset.psiPrev()(base_ip1_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[2], dx[axis_x],
140 : use_gradient_limiting);
141 4993 : const ControlPoints<T>& control_points_zx_y = GALS::INTERPOLATION::get_control_points(
142 9986 : levelset.psiPrev()(base_i_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_j_kp1)[0],
143 9986 : levelset.psiPrev()(base_ip1_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[0], dx[axis_x],
144 : use_gradient_limiting);
145 4993 : const ControlPoints<T>& control_points_yx_y = GALS::INTERPOLATION::get_control_points(
146 9986 : levelset.psiPrev()(base_i_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[0],
147 9986 : levelset.psiPrev()(base_ip1_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[0], dx[axis_x],
148 : use_gradient_limiting);
149 4993 : const ControlPoints<T>& control_points_yz_xy = GALS::INTERPOLATION::get_control_points(
150 9986 : levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[0], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[3],
151 9986 : levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[0], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[3],
152 9986 : dx[axis_z], use_gradient_limiting);
153 4993 : const ControlPoints<T>& control_points_yz_y = GALS::INTERPOLATION::get_control_points(
154 9986 : levelset.psiPrev()(base_i_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[1],
155 9986 : levelset.psiPrev()(base_i_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[1], dx[axis_z],
156 : use_gradient_limiting);
157 4993 : const ControlPoints<T>& control_points_yx_z = GALS::INTERPOLATION::get_control_points(
158 9986 : levelset.psiPrev()(base_i_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_k)[2],
159 9986 : levelset.psiPrev()(base_ip1_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[2], dx[axis_x],
160 : use_gradient_limiting);
161 4993 : const ControlPoints<T>& control_points_yzx_y = GALS::INTERPOLATION::get_control_points(
162 9986 : levelset.psiPrev()(base_i_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[0],
163 9986 : levelset.psiPrev()(base_ip1_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[0], dx[axis_x],
164 : use_gradient_limiting);
165 4993 : const ControlPoints<T>& control_points_yzx_z = GALS::INTERPOLATION::get_control_points(
166 9986 : levelset.psiPrev()(base_i_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_i_jp1_kp1)[2],
167 9986 : levelset.psiPrev()(base_ip1_jp1_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[2], dx[axis_x],
168 : use_gradient_limiting);
169 4993 : const ControlPoints<T>& control_points_xz_xy = GALS::INTERPOLATION::get_control_points(
170 9986 : levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[0], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[3],
171 9986 : levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[0], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[3],
172 9986 : dx[axis_z], use_gradient_limiting);
173 4993 : const ControlPoints<T>& control_points_xz_y = GALS::INTERPOLATION::get_control_points(
174 9986 : levelset.psiPrev()(base_ip1_j_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[1],
175 9986 : levelset.psiPrev()(base_ip1_j_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[1], dx[axis_z],
176 : use_gradient_limiting);
177 4993 : const ControlPoints<T>& control_points_xyz_xy = GALS::INTERPOLATION::get_control_points(
178 9986 : levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[0], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[3],
179 9986 : levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[0], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[3],
180 9986 : dx[axis_z], use_gradient_limiting);
181 4993 : const ControlPoints<T>& control_points_xyz_y = GALS::INTERPOLATION::get_control_points(
182 9986 : levelset.psiPrev()(base_ip1_jp1_k)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[1],
183 9986 : levelset.psiPrev()(base_ip1_jp1_kp1)[axis_y], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_kp1)[1], dx[axis_z],
184 : use_gradient_limiting);
185 4993 : const ControlPoints<T>& control_points_xy_z = GALS::INTERPOLATION::get_control_points(
186 9986 : levelset.psiPrev()(base_ip1_j_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_k)[1],
187 9986 : levelset.psiPrev()(base_ip1_jp1_k)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_jp1_k)[1], dx[axis_y],
188 : use_gradient_limiting);
189 4993 : const ControlPoints<T>& control_points_zxy_z = GALS::INTERPOLATION::get_control_points(
190 9986 : levelset.psiPrev()(base_ip1_j_kp1)[axis_z], levelset.phiMixedDerivativesPrev()(base_ip1_j_kp1)[1],
191 9986 : 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 4993 : const T& hx_by_three = dx[axis_x] * (T)one_third;
271 4993 : const T& hy_by_three = dx[axis_y] * (T)one_third;
272 4993 : const T& hz_by_three = dx[axis_z] * (T)one_third;
273 :
274 24965 : T bx[] = {B0(eta[0]), B1(eta[0]), B2(eta[0]), B3(eta[0])};
275 24965 : T by[] = {B0(eta[1]), B1(eta[1]), B2(eta[1]), B3(eta[1])};
276 24965 : T bz[] = {B0(eta[2]), B1(eta[2]), B2(eta[2]), B3(eta[2])};
277 24965 : T bx_prime[] = {B0_Prime(eta[0]), B1_Prime(eta[0]), B2_Prime(eta[0]), B3_Prime(eta[0])};
278 24965 : T by_prime[] = {B0_Prime(eta[1]), B1_Prime(eta[1]), B2_Prime(eta[1]), B3_Prime(eta[1])};
279 24965 : T bz_prime[] = {B0_Prime(eta[2]), B1_Prime(eta[2]), B2_Prime(eta[2]), B3_Prime(eta[2])};
280 :
281 164769 : 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 4993 : c_21_y + (dx[axis_z] * (T)one_third * c_21_y_z),
288 4993 : c_21_zy - (dx[axis_z] * (T)one_third * c_21_zy_z),
289 : c_21_zy,
290 : c_12_y,
291 4993 : c_12_y + (dx[axis_z] * (T)one_third * c_12_y_z),
292 4993 : 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 4993 : c_21_z + (dx[axis_x] * (T)one_third * c_21_z_x),
300 4993 : c_12_z + (dx[axis_x] * (T)one_third * c_12_z_x),
301 : c_21_zx,
302 4993 : c_21_x + (dx[axis_y] * (T)one_third * c_21_x_y),
303 4993 : 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 4993 : 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 4993 : c_21_zx + (hy_by_three * c_21_zx_y),
306 4993 : c_21_yx - (hy_by_three * c_21_yx_y),
307 4993 : 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 4993 : 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 4993 : c_21_yzx - (hy_by_three * c_21_yzx_y),
310 : c_21_yx,
311 4993 : c_21_yx + (hz_by_three * c_21_yx_z),
312 4993 : c_21_yzx - (hz_by_three * c_21_yzx_z),
313 : c_21_yzx,
314 : c_12_x,
315 4993 : c_12_x + (hz_by_three * c_12_x_z),
316 4993 : c_12_zx - (hz_by_three * c_12_zx_z),
317 : c_12_zx,
318 4993 : c_12_x + (hy_by_three * c_12_x_y),
319 4993 : 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 4993 : 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 4993 : c_12_zx + (hy_by_three * c_12_zx_y),
322 4993 : c_12_yx - (hy_by_three * c_12_yx_y),
323 4993 : 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 4993 : 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 4993 : c_12_yzx - (hy_by_three * c_12_yzx_y),
326 : c_12_yx,
327 4993 : c_12_yx + (hz_by_three * c_12_yx_z),
328 4993 : 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 4993 : c_21_xy + (hz_by_three * c_21_xy_z),
336 4993 : c_21_zxy - (hz_by_three * c_21_zxy_z),
337 : c_21_zxy,
338 : c_12_xy,
339 4993 : c_12_xy + (hz_by_three * c_12_xy_z),
340 4993 : 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 44937 : for (int i = 0, l = 0; i < 4; ++i)
348 179748 : for (int j = 0; j < 4; ++j)
349 718992 : for (int k = 0; k < 4; ++k) {
350 319552 : hermite_fields.phi_interpolated += bx[i] * by[j] * bz[k] * control_points_all[l];
351 319552 : ++l;
352 : }
353 :
354 44937 : for (int i = 0, l = 0; i < 4; ++i)
355 179748 : for (int j = 0; j < 4; ++j)
356 718992 : for (int k = 0; k < 4; ++k) {
357 319552 : hermite_fields.psi_interpolated[axis_x] += bx_prime[i] * by[j] * bz[k] * control_points_all[l];
358 319552 : ++l;
359 : }
360 4993 : hermite_fields.psi_interpolated[axis_x] *= one_over_dx[axis_x];
361 :
362 44937 : for (int i = 0, l = 0; i < 4; ++i)
363 179748 : for (int j = 0; j < 4; ++j)
364 718992 : for (int k = 0; k < 4; ++k) {
365 319552 : hermite_fields.psi_interpolated[axis_y] += bx[i] * by_prime[j] * bz[k] * control_points_all[l];
366 319552 : ++l;
367 : }
368 4993 : hermite_fields.psi_interpolated[axis_y] *= one_over_dx[axis_y];
369 :
370 44937 : for (int i = 0, l = 0; i < 4; ++i)
371 179748 : for (int j = 0; j < 4; ++j)
372 718992 : for (int k = 0; k < 4; ++k) {
373 319552 : hermite_fields.psi_interpolated[axis_z] += bx[i] * by[j] * bz_prime[k] * control_points_all[l];
374 319552 : ++l;
375 : }
376 4993 : hermite_fields.psi_interpolated[axis_z] *= one_over_dx[axis_z];
377 :
378 4993 : return hermite_fields;
379 : }
380 :
381 : template <typename T>
382 4 : 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 8 : const GALS::CPU::Vec3<int> num_cells_interp = x_interp.numCells();
389 4 : const T_GRID& grid = levelset.grid();
390 8 : 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 90 : for (int i = 0; i < num_cells_interp[0]; ++i)
394 969 : for (int j = 0; j < num_cells_interp[1]; ++j)
395 10449 : for (int k = 0; k < num_cells_interp[2]; ++k) {
396 4993 : const auto& hermite_fields = this->interpolate(grid, x_interp(i, j, k), levelset);
397 :
398 4993 : levelset.phiInterpPrev()(i, j, k) = hermite_fields.phi_interpolated;
399 4993 : levelset.psiInterpPrev()(i, j, k) = hermite_fields.psi_interpolated;
400 : }
401 4 : }
402 :
403 3 : template class GALS::INTERPOLATION::Hermite<double, GALS::CPU::Grid<double, 3>>;
|
