OpenLB 1.7
Loading...
Searching...
No Matches
blockStructure.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2021 Mathias Krause, Adrian Kummerlaender
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 BLOCK_STRUCTURE_H
25#define BLOCK_STRUCTURE_H
26
27#include <cstdint>
28#include <type_traits>
29
30#include "vector.h"
31#include "olbDebug.h"
32
33namespace olb {
34
36using CellID = std::uint32_t;
37
39using CellDistance = std::int64_t;
40
42template <unsigned D>
44
46template <typename T, unsigned D>
48
50
53template <unsigned D>
55protected:
59
61
62public:
63 static_assert(D == 2 || D == 3, "Only D=2 and D=3 are supported");
64
65 BlockStructureD(Vector<int,D> size, int padding=0):
66 _core(size),
67 _size(size + 2*padding),
68 _padding(padding)
69 {
70 if constexpr (D == 3) {
71 _projection = {_size[1]*_size[2], _size[2], 1};
72 } else {
73 _projection = {_size[1], 1};
74 }
75
76 if (getNcells()-1 > std::numeric_limits<CellID>::max()) {
77 throw std::invalid_argument("Cell count must not exceed cell index space");
78 }
79 };
80
82 BlockStructureD(1, 0)
83 { };
84
86 int getNx() const
87 {
88 return _core[0];
89 };
91 int getNy() const
92 {
93 return _core[1];
94 };
96 int getNz() const
97 {
98 static_assert(D == 3, "z-component only available in 3D");
99 return _core[2];
100 };
101
103 {
104 return _core;
105 }
106
108 int getPadding() const
109 {
110 return _padding;
111 };
112
114 std::size_t getNcells() const
115 {
116 if constexpr (D == 3) {
117 return _size[0] * _size[1] * _size[2];
118 } else {
119 return _size[0] * _size[1];
120 }
121 __builtin_unreachable();
122 }
123
126 {
127 OLB_PRECONDITION(isInside(latticeR));
128 return (latticeR+_padding) * _projection;
129 }
130
131 template <typename... L>
132 std::enable_if_t<sizeof...(L) == D, CellID>
133 getCellId(L... latticeR) const
134 {
135 return LatticeR<D>{latticeR+_padding...} * _projection;
136 }
137
140 {
141 return dir * _projection;
142 }
143
145 bool isInside(LatticeR<D> latticeR) const
146 {
147 return latticeR >= -_padding && latticeR < _core + _padding;
148 };
149
151 bool isInsideCore(LatticeR<D> latticeR) const
152 {
153 return latticeR >= 0 && latticeR < _core;
154 };
155
157 bool isPadding(LatticeR<D> latticeR) const
158 {
159 return isInside(latticeR) && !(latticeR >= 0 && latticeR < _core);
160 };
161
162 template <typename... L>
163 std::enable_if_t<sizeof...(L) == D, bool>
164 isInside(L... latticeR) const
165 {
166 return isInside({latticeR...});
167 };
168
171 {
172 auto lower = latticeR + _padding;
173 auto upper = LatticeR<D>([&](unsigned iDim) -> int {
174 int x = lower[iDim] - _size[iDim] + 1;
175 return x < 0 ? -x : 0;
176 });
177 auto y = std::min(*std::min_element(lower.data(), lower.data() + D),
178 *std::min_element(upper.data(), upper.data() + D));
179 return y;
180 };
181
182 template <typename F>
183 void forSpatialLocations(F f) const
184 {
185 using loc = typename LatticeR<D>::value_t;
186 for (loc iX=-_padding; iX < _core[0] + _padding; ++iX) {
187 for (loc iY=-_padding; iY < _core[1] + _padding; ++iY) {
188 if constexpr (D == 3) {
189 for (loc iZ=-_padding; iZ < _core[2] + _padding; ++iZ) {
190 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
191 f({iX,iY,iZ});
192 } else {
193 f(iX,iY,iZ);
194 }
195 }
196 } else {
197 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
198 f({iX,iY});
199 } else {
200 f(iX,iY);
201 }
202 }
203 }
204 }
205 };
206
207 template <typename F>
209 {
210 using loc = typename LatticeR<D>::value_t;
211 #ifdef PARALLEL_MODE_OMP
212 #pragma omp parallel for schedule(dynamic,1)
213 #endif
214 for (loc iX=-_padding; iX < _core[0] + _padding; ++iX) {
215 for (loc iY=-_padding; iY < _core[1] + _padding; ++iY) {
216 if constexpr (D == 3) {
217 for (loc iZ=-_padding; iZ < _core[2] + _padding; ++iZ) {
218 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
219 f({iX,iY,iZ});
220 } else {
221 f(iX,iY,iZ);
222 }
223 }
224 } else {
225 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
226 f({iX,iY});
227 } else {
228 f(iX,iY);
229 }
230 }
231 }
232 }
233 }
234
235 template <typename F>
237 {
238 using loc = typename LatticeR<D>::value_t;
239 for (loc iX=std::max(-_padding, min[0]); iX < std::min(_core[0] + _padding, max[0]+1); ++iX) {
240 for (loc iY=std::max(-_padding, min[1]); iY < std::min(_core[1] + _padding, max[1]+1); ++iY) {
241 if constexpr (D == 3) {
242 for (loc iZ=std::max(-_padding, min[2]); iZ < std::min(_core[2] + _padding, max[2]+1); ++iZ) {
243 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
244 f({iX,iY,iZ});
245 } else {
246 f(iX,iY,iZ);
247 }
248 }
249 } else {
250 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
251 f({iX,iY});
252 } else {
253 f(iX,iY);
254 }
255 }
256 }
257 }
258 };
259
260 template <typename F>
262 {
263 using loc = typename LatticeR<D>::value_t;
264 for (loc iX=0; iX < _core[0]; ++iX) {
265 for (loc iY=0; iY < _core[1]; ++iY) {
266 if constexpr (D == 3) {
267 for (loc iZ=0; iZ < _core[2]; ++iZ) {
268 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
269 f({iX,iY,iZ});
270 } else {
271 f(iX,iY,iZ);
272 }
273 }
274 } else {
275 if constexpr (std::is_invocable_v<F, LatticeR<D>>) {
276 f({iX,iY});
277 } else {
278 f(iX,iY);
279 }
280 }
281 }
282 }
283 };
284
285 template <typename F>
286 void forCellIndices(F f) const
287 {
288 for (CellID iCell=0; iCell < getNcells(); ++iCell) {
289 f(iCell);
290 }
291 }
292
293};
294
295template <typename DESCRIPTOR>
297
298
299}
300
301#endif
Base of a regular block.
CellDistance getNeighborDistance(LatticeR< D > dir) const
Get 1D neighbor distance.
void forCellIndices(F f) const
std::size_t getNcells() const
Get number of cells.
bool isInsideCore(LatticeR< D > latticeR) const
Return whether location is inside core.
BlockStructureD(Vector< int, D > size, int padding=0)
int getNy() const
Read only access to block height.
void forSpatialLocations(F f) const
void forSpatialLocationsParallel(F f) const
CellDistance getNeighborhoodRadius(LatticeR< D > latticeR) const
Return maximum valid neighborhood sphere radius w.r.t. latticeR.
int getPadding() const
Read only access to padding.
void forCoreSpatialLocations(F f) const
LatticeR< D > _projection
bool isPadding(LatticeR< D > latticeR) const
Return whether location is valid.
int getNx() const
Read only access to block width.
CellID getCellId(LatticeR< D > latticeR) const
Get 1D cell ID.
int getNz() const
Read only access to block height.
void forSpatialLocations(LatticeR< D > min, LatticeR< D > max, F f) const
bool isInside(LatticeR< D > latticeR) const
Return whether location is valid.
std::enable_if_t< sizeof...(L)==D, bool > isInside(L... latticeR) const
std::enable_if_t< sizeof...(L)==D, CellID > getCellId(L... latticeR) const
LatticeR< D > getExtent() const
Plain old scalar vector.
Definition vector.h:47
Top level namespace for all of OpenLB.
std::uint32_t CellID
Type for sequential block-local cell indices.
std::int64_t CellDistance
Type for in-memory distance of block-local cell indices.
#define OLB_PRECONDITION(COND)
Definition olbDebug.h:46
efficient implementation of a vector class