OpenLB 1.8.1
Loading...
Searching...
No Matches
interpolation2d.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2023-24 Julius Jessberger
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24#ifndef INTERPOLATION_2D_H
25#define INTERPOLATION_2D_H
26
27#include "core/vector.h"
28
29namespace olb {
30
36// Caution: Developer has to guarantee that the data are accessible.
37// It is expected that the point of interest (=cell+distance) lies nearby cell
38template <unsigned DataDim, typename CELL, typename RESULT, typename F, typename V = typename CELL::value_t>
39void interpolate2d(CELL& cell, RESULT& result, const Vector<V,2> distance, F f) any_platform
40{
41 const Vector<int,2> floorV (util::floor(distance[0]), util::floor(distance[1]));
42 Vector<Vector<int,2>,4> surroundingPoints (floorV);
43 surroundingPoints[1][0] += 1;
44 surroundingPoints[2][1] += 1;
45 surroundingPoints[3][0] += 1; surroundingPoints[3][1] += 1;
46
47 V resNeighbor[DataDim];
48 for (unsigned i=0; i<DataDim; ++i) {
49 result[i] = V();
50 }
51 for (auto point : surroundingPoints) {
52 const Vector<V,2> dist = distance - point;
53 const V weight = (V(1) - util::abs(dist[0])) * (V(1) - util::abs(dist[1]));
54 f(cell.neighbor(point), resNeighbor);
55 for (unsigned i = 0; i < DataDim; ++i) {
56 result[i] += weight * resNeighbor[i];
57 }
58 }
59}
60
62template <unsigned DataDim, typename CELL, typename RESULT, typename F, typename V = typename CELL::value_t>
63void interpolate2d_help(CELL& cell, RESULT& result, const Vector<V,2> distance, F f,
64 Vector<int,2> pointA, Vector<int,2> pointB, Vector<int,2> pointC) any_platform
65{
66 // express point as a linear combination of the two sides
67 const auto lc = util::solveLinearSystem(Vector<V,2>(pointB - pointA),
68 Vector<V,2>(pointC - pointA),
69 distance - pointA);
70 // evaluate linear function
71 V resNeighbor[DataDim];
72 for (unsigned i=0; i<DataDim; ++i) {
73 result[i] = V();
74 }
75 f(cell.neighbor(pointA), resNeighbor);
76 for (unsigned i = 0; i < DataDim; ++i) {
77 result[i] += (1 - lc[0] - lc[1]) * resNeighbor[i];
78 }
79 f(cell.neighbor(pointB), resNeighbor);
80 for (unsigned i = 0; i < DataDim; ++i) {
81 result[i] += lc[0] * resNeighbor[i];
82 }
83 f(cell.neighbor(pointC), resNeighbor);
84 for (unsigned i = 0; i < DataDim; ++i) {
85 result[i] += lc[1] * resNeighbor[i];
86 }
87}
88
90template <unsigned DataDim, typename CELL, typename RESULT, typename F, typename V = typename CELL::value_t>
91void interpolate2d_help(CELL& cell, RESULT& result, const Vector<V,2> distance, F f,
92 Vector<int,2> pointA, Vector<int,2> pointB) any_platform
93{
94 const V refDistanceSq = norm_squared(Vector<V,2>(pointA - pointB));
95 const V alpha = ((distance - pointB) * (pointA - pointB)) / refDistanceSq;
96 V resNeighbor[DataDim];
97 for (unsigned i=0; i<DataDim; ++i) {
98 result[i] = V();
99 }
100 f(cell.neighbor(pointA), resNeighbor);
101 for (unsigned i = 0; i < DataDim; ++i) {
102 result[i] += alpha * resNeighbor[i];
103 }
104 f(cell.neighbor(pointB), resNeighbor);
105 for (unsigned i = 0; i < DataDim; ++i) {
106 result[i] += (1 - alpha) * resNeighbor[i];
107 }
108}
109
119// Caution: Developer has to guarantee that the data are accessible.
120// It is expected that the point of interest (=cell+distance) lies nearby cell
121template <unsigned NEIGHBORS, unsigned DataDim, typename CELL, typename RESULT, typename F, typename V = typename CELL::value_t>
122void interpolate2d(CELL& cell, RESULT& result, const Vector<V,2> distance, F f) any_platform
123{
124 // four neighbors: bilinear interpolation
125 if constexpr (NEIGHBORS == 1111) {
126 interpolate2d<DataDim,CELL,RESULT,F,V>(cell, result, distance, f);
127 }
128 // three neighbors: linear interpolation
129 else if constexpr (NEIGHBORS == 1110) {
130 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1]));
131 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1]));
132 const Vector<int,2> pointC (util::floor(distance[0]), util::floor(distance[1])+1);
133 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB, pointC);
134 } else if constexpr (NEIGHBORS == 1101) {
135 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1]));
136 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1]));
137 const Vector<int,2> pointC (util::floor(distance[0]+1), util::floor(distance[1])+1);
138 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB, pointC);
139 } else if constexpr (NEIGHBORS == 1011) {
140 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1]));
141 const Vector<int,2> pointB (util::floor(distance[0]), util::floor(distance[1])+1);
142 const Vector<int,2> pointC (util::floor(distance[0]+1), util::floor(distance[1])+1);
143 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB, pointC);
144 } else if constexpr (NEIGHBORS == 111) {
145 const Vector<int,2> pointA (util::floor(distance[0]+1), util::floor(distance[1]));
146 const Vector<int,2> pointB (util::floor(distance[0]), util::floor(distance[1])+1);
147 const Vector<int,2> pointC (util::floor(distance[0]+1), util::floor(distance[1])+1);
148 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB, pointC);
149 }
150 else if constexpr (NEIGHBORS == 1100) {
151 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1]));
152 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1]));
153 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB);
154 } else if constexpr (NEIGHBORS == 11) {
155 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1])+1);
156 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1])+1);
157 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB);
158 } else if constexpr (NEIGHBORS == 1010) {
159 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1]));
160 const Vector<int,2> pointB (util::floor(distance[0]), util::floor(distance[1])+1);
161 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB);
162 } else if constexpr (NEIGHBORS == 101) {
163 const Vector<int,2> pointA (util::floor(distance[0]+1), util::floor(distance[1]));
164 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1])+1);
165 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB);
166 } else if constexpr (NEIGHBORS == 1001) {
167 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1]));
168 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1])+1);
169 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB);
170 } else if constexpr (NEIGHBORS == 110) {
171 const Vector<int,2> pointA (util::floor(distance[0]), util::floor(distance[1])+1);
172 const Vector<int,2> pointB (util::floor(distance[0]+1), util::floor(distance[1]));
173 interpolate2d_help<DataDim,CELL,RESULT,F,V>(cell, result, distance, f, pointA, pointB);
174 }
175 // one neighbor: constant extrapolation
176 else if constexpr (NEIGHBORS == 1) {
177 const Vector<int,2> point (util::floor(distance[0])+1, util::floor(distance[1])+1);
178 f(cell.neighbor(point), result);
179 } else if constexpr (NEIGHBORS == 10) {
180 const Vector<int,2> point (util::floor(distance[0]), util::floor(distance[1])+1);
181 f(cell.neighbor(point), result);
182 } else if constexpr (NEIGHBORS == 100) {
183 const Vector<int,2> point (util::floor(distance[0])+1, util::floor(distance[1]));
184 f(cell.neighbor(point), result);
185 } else if constexpr (NEIGHBORS == 1000) {
186 const Vector<int,2> point (util::floor(distance[0]), util::floor(distance[1]));
187 f(cell.neighbor(point), result);
188 }
189 else {
190 throw std::invalid_argument("Invalid NEIGHBORS");
191 }
192}
193
201// Caution: Developer has to guarantee that the data are accessible.
202// It is expected that the point of interest lies nearby cell
203template <unsigned NEIGHBORS, typename FIELD, typename CELL, typename RESULT, typename V = typename CELL::value_t>
204void interpolate2d(CELL& cell, RESULT& result, const Vector<V,2> distance) any_platform
205{
206 FieldD<V, typename CELL::descriptor_t, FIELD> v;
207 constexpr unsigned DataDim = CELL::descriptor_t::template size<FIELD>();
208 auto f = [DataDim](const auto& c, RESULT& res) -> void {
209 for (unsigned i = 0; i < DataDim; ++i) {
210 res[i] = c.template getFieldComponent<FIELD>(i);
211 }
212 };
213 interpolate2d<NEIGHBORS,DataDim,CELL,RESULT,decltype(f),V>(cell, result, distance, f);
214}
215
219// Caution: Developer has to guarantee that the data are accessible.
220// It is expected that the point of interest lies nearby cell
221template <typename FIELD, typename CELL, typename RESULT, typename V = typename CELL::value_t>
222void interpolate2d(CELL& cell, RESULT& result, const Vector<V,2> distance) any_platform
223{
224 interpolate2d<1111,FIELD,CELL,RESULT,V>(cell, result, distance);
225}
226
227
228}
229
230#endif
olb::Vector
Plain old scalar vector.
Definition vectorHelpers.h:41
olb::util::abs
ADf< T, DIM > abs(const ADf< T, DIM > &a)
Definition aDiff.h:1019
olb::util::floor
ADf< T, DIM > floor(const ADf< T, DIM > &a)
Definition aDiff.h:869
olb::util::solveLinearSystem
auto solveLinearSystem(const Matrix< T, 2, 2 > &a, const Vector< T, 2 > &rhs)
Solve a * x = rhs.
Definition matrix.h:158
olb
Top level namespace for all of OpenLB.
Definition advectionDiffusionDirichlet.h:27
olb::norm_squared
constexpr T norm_squared(const ScalarVector< T, D, IMPL > &a) any_platform
Squared euclidean vector norm.
Definition scalarVector.h:55
olb::interpolate2d
void interpolate2d(CELL &cell, RESULT &result, const Vector< V, 2 > distance, F f) any_platform
Bilinear interpolation between four neighboring mesh points.
Definition interpolation2d.h:39
olb::interpolate2d_help
void interpolate2d_help(CELL &cell, RESULT &result, const Vector< V, 2 > distance, F f, Vector< int, 2 > pointA, Vector< int, 2 > pointB, Vector< int, 2 > pointC) any_platform
Linear interpolation between three points.
Definition interpolation2d.h:63
any_platform
#define any_platform
Define preprocessor macros for device-side functions, constant storage.
Definition platform.h:77
vector.h
efficient implementation of a vector class
Generated on Sat May 17 2025 10:48:04 for OpenLB by doxygen 1.12.0