OpenLB 1.7
Loading...
Searching...
No Matches
cuboid2D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2007, 2014 Mathias J. Krause
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
29#ifndef CUBOID_2D_HH
30#define CUBOID_2D_HH
31
32#include <iostream>
33#include <math.h>
34#include <vector>
35#include "cuboid2D.h"
36#include "dynamics/lbm.h"
37#include "io/ostreamManager.h"
38
39
40namespace olb {
41
43
44template<typename T>
45Cuboid2D<T>::Cuboid2D(T globPosX, T globPosY, T delta,int nX, int nY)
46 : _weight(std::numeric_limits<size_t>::max()), clout(std::cout,"Cuboid2D")
47{
48 init(globPosX, globPosY, delta, nX, nY);
49}
50
51template<typename T>
53 : _weight(std::numeric_limits<size_t>::max()), clout(std::cout,"Cuboid2D")
54{
55 this->init(origin[0], origin[1], delta, extend[0], extend[1]);
56}
57
58// copy constructor
59template<typename T>
60Cuboid2D<T>::Cuboid2D(Cuboid2D<T> const& rhs, int overlap) : clout(std::cout,"Cuboid2D")
61{
62 this->init(rhs._globPosX - rhs._delta*overlap, rhs._globPosY - rhs._delta*overlap,
63 rhs._delta, rhs._nX + 2*overlap, rhs._nY + 2*overlap);
64 _weight = rhs._weight;
65}
66
67// copy assignment
68template<typename T>
70{
71 this->init(rhs._globPosX, rhs._globPosY, rhs._delta, rhs._nX, rhs._nY);
72 _weight = rhs._weight;
73 return *this;
74}
75
76template<typename T>
77void Cuboid2D<T>::init(T globPosX, T globPosY, T delta, int nX, int nY)
78{
79 _globPosX = globPosX;
80 _globPosY = globPosY;
81 _delta = delta;
82 _nX = nX;
83 _nY = nY;
84}
85
86
87template<typename T>
89{
90 return _globPosX;
91}
92
93template<typename T>
95{
96 return _globPosY;
97}
98
99template<typename T>
101{
102 return Vector<T,2> (_globPosX,_globPosY);
103}
104
105template<typename T>
107{
108 return _delta;
109}
110
111template<typename T>
113{
114 return _nX;
115}
116
117template<typename T>
119{
120 return _nY;
121}
122
123template<typename T>
125{
126 return Vector<int,2> (_nX, _nY);
127}
128
129template<typename T>
131{
132 return _nY*_nX*_delta*_delta;
133}
134
135template<typename T>
137{
138 return static_cast<size_t>(_nY)*static_cast<size_t>(_nX);
139}
140
141template<typename T>
143{
144 return 2*_nY*_delta + 2*_nX*_delta;
145}
146
147template<typename T>
149{
150 return 2*_nY + 2*_nX -4;
151}
152
153template<typename T>
155{
156 clout << "--------Cuboid Details----------" << std::endl;
157 clout << " Left Corner (x/y): " << "\t" << "(" << this->get_globPosX() << "/" << this->get_globPosY() << ")" << std::endl;
158 clout << " Delta: " << "\t" << "\t" << this->getDeltaR() << std::endl;
159 clout << " Perimeter: " << "\t" << "\t" << this->getPhysPerimeter() << std::endl;
160 clout << " Volume: " << "\t" << "\t" << this->getPhysVolume() << std::endl;
161 clout << " Extent (x/y): " << "\t" << "\t" << "(" << this->getNx() << "/" << this->getNy() << ")" << std::endl;
162 clout << " Nodes at Perimeter: " << "\t" << this->getLatticePerimeter() << std::endl;
163 clout << " Nodes in Volume: " << "\t" << this->getLatticeVolume() << std::endl;
164 clout << "--------------------------------" << std::endl;
165}
166
167template<typename T>
168void Cuboid2D<T>::getPhysR(T physR[2], const int latticeR[2]) const
169{
170 physR[0] = _globPosX + latticeR[0]*_delta;
171 physR[1] = _globPosY + latticeR[1]*_delta;
172}
173
174template<typename T>
175void Cuboid2D<T>::getPhysR(T physR[2], const int& iX, const int& iY) const
176{
177 physR[0] = _globPosX + iX*_delta;
178 physR[1] = _globPosY + iY*_delta;
179}
180
181template<typename T>
182void Cuboid2D<T>::getPhysR(T physR[2], LatticeR<2> latticeR) const
183{
184 physR[0] = _globPosX + latticeR[0]*_delta;
185 physR[1] = _globPosY + latticeR[1]*_delta;
186}
187
188template<typename T>
189void Cuboid2D<T>::getLatticeR(int latticeR[2], const T physR[2]) const
190{
191 latticeR[0] = (int)util::floor( (physR[0] - _globPosX )/_delta +.5);
192 latticeR[1] = (int)util::floor( (physR[1] - _globPosY )/_delta +.5);
193}
194
195template<typename T>
196void Cuboid2D<T>::getLatticeR(int latticeR[2], const Vector<T,2>& physR) const
197{
198 latticeR[0] = (int)util::floor( (physR[0] - _globPosX )/_delta +.5);
199 latticeR[1] = (int)util::floor( (physR[1] - _globPosY )/_delta +.5);
200}
201
202template<typename T>
203bool Cuboid2D<T>::checkPoint(T globX, T globY, int overlap) const
204{
205
206 if (_globPosX <= globX + T(0.5 + overlap)*_delta &&
207 _globPosX + T(_nX+overlap)*_delta > globX + 0.5*_delta &&
208 _globPosY <= globY + T(0.5 + overlap)*_delta &&
209 _globPosY + T(_nY+overlap)*_delta > globY + 0.5*_delta ) {
210 return true;
211 }
212 else {
213 return false;
214 }
215}
216
217template<typename T>
218bool Cuboid2D<T>::checkPoint(Vector<T,2>& globXY, int overlap) const
219{
220 return checkPoint(globXY[0], globXY[1], overlap);
221}
222
223template<typename T>
224bool Cuboid2D<T>::checkPoint(T globX, T globY, int &locX, int &locY, int overlap) const
225{
226 if (overlap!=0) {
227 Cuboid2D tmp(_globPosX - overlap*_delta, _globPosY - overlap*_delta, _delta, _nX + overlap*2, _nY + overlap*2);
228 return tmp.checkPoint(globX, globY, locX, locY);
229 }
230 else if (!checkPoint(globX, globY)) {
231 return false;
232 }
233 else {
234 locX = (int)util::floor((globX - (T)_globPosX)/_delta + .5);
235 locY = (int)util::floor((globY - (T)_globPosY)/_delta + .5);
236 return true;
237 }
238}
239
240template<typename T>
241bool Cuboid2D<T>::checkInters(T globX0, T globX1, T globY0, T globY1, int overlap) const
242{
243
244 T locX0d = util::max(_globPosX-overlap*_delta,globX0);
245 T locY0d = util::max(_globPosY-overlap*_delta,globY0);
246 T locX1d = util::min(_globPosX+(_nX+overlap-1)*_delta,globX1);
247 T locY1d = util::min(_globPosY+(_nY+overlap-1)*_delta,globY1);
248
249 if (!(locX1d>=locX0d && locY1d>=locY0d)) {
250 return false;
251 }
252 return true;
253}
254
255template<typename T>
256bool Cuboid2D<T>::checkInters(T globX, T globY, int overlap) const
257{
258 return checkInters(globX, globX, globY, globY, overlap);
259}
260
261template<typename T>
262bool Cuboid2D<T>::checkInters(T globX0, T globX1, T globY0, T globY1,
263 int &locX0, int &locX1, int &locY0, int &locY1, int overlap) const
264{
265 if (overlap!=0) {
266 Cuboid2D tmp(_globPosX - overlap*_delta, _globPosY - overlap*_delta, _delta, _nX + overlap*2, _nY + overlap*2);
267 return tmp.checkInters(globX0, globX1, globY0, globY1, locX0, locX1, locY0, locY1);
268 }
269 else if (!checkInters(globX0, globX1, globY0, globY1)) {
270 locX0 = 1;
271 locX1 = 0;
272 locY0 = 1;
273 locY1 = 0;
274 return false;
275 }
276 else {
277 locX0 = 0;
278 for (int i=0; _globPosX + i*_delta < globX0; i++) {
279 locX0 = i+1;
280 }
281 locX1 = _nX-1;
282 for (int i=_nX-1; _globPosX + i*_delta > globX1; i--) {
283 locX1 = i-1;
284 }
285 locY0 = 0;
286 for (int i=0; _globPosY + i*_delta < globY0; i++) {
287 locY0 = i+1;
288 }
289 locY1 = _nY-1;
290 for (int i=_nY-1; _globPosY + i*_delta > globY1; i--) {
291 locY1 = i-1;
292 }
293 return true;
294 }
295}
296
297template<typename T>
298void Cuboid2D<T>::divide(int nX, int nY, std::vector<Cuboid2D<T> > &childrenC) const
299{
300 T globPosX_child, globPosY_child;
301 int xN_child = 0;
302 int yN_child = 0;
303
304 globPosX_child = _globPosX;
305 globPosY_child = _globPosY;
306
307 for (int iX=0; iX<nX; iX++) {
308 for (int iY=0; iY<nY; iY++) {
309 xN_child = (_nX+nX-iX-1)/nX;
310 yN_child = (_nY+nY-iY-1)/nY;
311 Cuboid2D<T> child(globPosX_child, globPosY_child, _delta, xN_child, yN_child);
312 childrenC.push_back(child);
313 globPosY_child += yN_child*_delta;
314 }
315 globPosY_child = _globPosY;
316 globPosX_child += xN_child*_delta;
317 }
318}
319
320template<typename T>
321void Cuboid2D<T>::divide(int p, std::vector<Cuboid2D<T> > &childrenC) const
322{
323
324 OLB_PRECONDITION(p>0);
325 int nX = 0;
326 int nY = 0;
327 T ratio;
328 T bestRatio = (T)_nX/(T)_nY;
329 T difRatio = util::fabs(bestRatio - 1) + 1;
330 for (int i=1; i<= p; i++) {
331 int j = p / i;
332 if (i*j<=p) {
333 if ( util::fabs(bestRatio - (T)i/(T)j) <= difRatio) {
334 difRatio = util::fabs(bestRatio - (T)i/(T)j);
335 nX = i;
336 nY = j;
337 }
338 }
339 }
340
341 ratio = T(nX)/(T)nY;
342 int rest = p - nX*nY;
343
344 if (rest==0) {
345 divide(nX,nY,childrenC);
346 return;
347 }
348
349 if (ratio < bestRatio && (nY-rest) >= 0) {
350 int n_QNoInsertions = nX*(nY-rest);
351 T bestVolume_QNoInsertions = (T)_nX*_nY * n_QNoInsertions/(T)p;
352 int yN_QNoInsertions = (int)(bestVolume_QNoInsertions / (T)_nX);
353 int xN_QNoInsertions = _nX;
354 int yN_QInsertions = _nY-yN_QNoInsertions;
355 int xN_QInsertions = _nX;
356 Cuboid2D<T> firstChildQ(_globPosX, _globPosY, _delta, xN_QNoInsertions, yN_QNoInsertions);
357 Cuboid2D<T> secondChildQ(_globPosX, _globPosY+yN_QNoInsertions*_delta, _delta, xN_QInsertions, yN_QInsertions);
358 firstChildQ.divide(nX, nY-rest, childrenC);
359 secondChildQ.divide(nX+1,rest, childrenC);
360 }
361 else {
362 int n_QNoInsertions = nY*(nX-rest);
363 T bestVolume_QNoInsertions = (T)_nX*_nY * n_QNoInsertions/(T)p;
364 int xN_QNoInsertions = (int)(bestVolume_QNoInsertions / (T)_nY + 0.9999);
365 int yN_QNoInsertions = _nY;
366 int xN_QInsertions = _nX-xN_QNoInsertions;
367 int yN_QInsertions = _nY;
368 Cuboid2D<T> firstChildQ(_globPosX, _globPosY, _delta, xN_QNoInsertions, yN_QNoInsertions);
369 Cuboid2D<T> secondChildQ(_globPosX+xN_QNoInsertions*_delta, _globPosY, _delta, xN_QInsertions, yN_QInsertions);
370 firstChildQ.divide(nX-rest, nY, childrenC);
371 secondChildQ.divide(rest,nY+1, childrenC);
372 }
373}
374
375
376template<typename T>
377void Cuboid2D<T>::resize(int iX, int iY, int nX, int nY)
378{
379 _globPosX = _globPosX+iX*_delta;
380 _globPosY = _globPosY+iY*_delta;
381 _nX = nX;
382 _nY = nY;
383}
384
385template<typename T>
387{
388 if (_weight == std::numeric_limits<size_t>::max()) {
389 return getLatticeVolume();
390 }
391 else {
392 return _weight;
393 }
394}
395
396template<typename T>
397void Cuboid2D<T>::setWeight(size_t fullCells)
398{
399 _weight = fullCells;
400}
401
402
403} // namespace olb
404
405#endif
A regular single 2D cuboid is the basic component of a 2D cuboid structure which defines the grid.
Definition cuboid2D.h:54
size_t getLatticeVolume() const
Returns the number of Nodes in the volume.
Definition cuboid2D.hh:136
void print() const
Prints cuboid details.
Definition cuboid2D.hh:154
T getPhysVolume() const
Returns the volume of the cuboid.
Definition cuboid2D.hh:130
int getLatticePerimeter() const
Returns the number of Nodes at the perimeter.
Definition cuboid2D.hh:148
Cuboid2D(T globPosX, T globPosY, T delta, int nX, int nY)
Construction of a cuboid.
Definition cuboid2D.hh:45
void setWeight(size_t fullCells)
Sets the number of full cells.
Definition cuboid2D.hh:397
void getLatticeR(int latticeR[2], const T physR[2]) const
Definition cuboid2D.hh:189
bool checkInters(T globX0, T globX1, T globY0, T globY1, int overlap=0) const
Checks whether there is an intersection with the cuboid extended with an layer of size overlap*delta.
Definition cuboid2D.hh:241
T get_globPosY() const
Definition cuboid2D.hh:94
size_t getWeight() const
Returns the number of full cells.
Definition cuboid2D.hh:386
void divide(int p, int q, std::vector< Cuboid2D< T > > &childrenC) const
Divides the cuboid in p*q cuboids and adds them to the given vector.
Definition cuboid2D.hh:298
Vector< int, 2 > const getExtent() const
Read only access to the number of voxels in every dimension.
Definition cuboid2D.hh:124
T getDeltaR() const
Read access to the distance of cuboid nodes.
Definition cuboid2D.hh:106
bool checkPoint(T globX, T globY, int overlap=0) const
Checks whether a point (globX/globY) is contained in the cuboid extended with an layer of size overla...
Definition cuboid2D.hh:203
void resize(int X, int Y, int nX, int nY)
resize the cuboid to the passed size
Definition cuboid2D.hh:377
void getPhysR(T physR[2], const int latticeR[2]) const
Definition cuboid2D.hh:168
Vector< T, 2 > const getOrigin() const
Read only access to left lower corner coordinates.
Definition cuboid2D.hh:100
int getNx() const
Read access to cuboid width.
Definition cuboid2D.hh:112
int getNy() const
Read access to cuboid height.
Definition cuboid2D.hh:118
T get_globPosX() const
Read access to left lower corner coordinates.
Definition cuboid2D.hh:88
T getPhysPerimeter() const
Returns the perimeter of the cuboid.
Definition cuboid2D.hh:142
Cuboid2D & operator=(Cuboid2D const &rhs)
Copy assignment.
Definition cuboid2D.hh:69
void init(T globPosX, T globPosY, T delta, int nX, int nY)
Initializes the cuboid.
Definition cuboid2D.hh:77
Plain old scalar vector.
Definition vector.h:47
The description of a single 2D cuboid – header file.
ADf< T, DIM > floor(const ADf< T, DIM > &a)
Definition aDiff.h:869
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
Top level namespace for all of OpenLB.
#define OLB_PRECONDITION(COND)
Definition olbDebug.h:46