OpenLB 1.7
Loading...
Searching...
No Matches
blockGeometry.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2011, 2014 Mathias J. Krause, Simon Zimny
4 * 2021 Clara Schragmann, Adrian Kummerlaender
5 * E-mail contact: info@openlb.net
6 * The most recent release of OpenLB can be downloaded at
7 * <http://www.openlb.net/>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public
20 * License along with this program; if not, write to the Free
21 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
23 */
24
29#ifndef BLOCK_GEOMETRY_HH
30#define BLOCK_GEOMETRY_HH
31
33
36
39
40namespace olb {
41
42
43template<typename T, unsigned D>
44BlockGeometry<T,D>::BlockGeometry(Cuboid<T,D>& cuboid, int padding, int iCglob)
45 : BlockStructureD<D>(cuboid.getExtent(), padding),
46 _data(this->getNcells()),
47 _communicatable(_data),
48 _cuboid(cuboid),
49 _iCglob(iCglob),
50 _statistics(this),
51 clout(std::cout, ("BlockGeometry" + std::to_string(D) + "D"))
52{
53 _statistics.update(false);
54}
55
56template<typename T, unsigned D>
58{
59 return _statistics;
60}
61
62template<typename T, unsigned D>
67
68template<typename T, unsigned D>
70{
71 return _cuboid.getOrigin();
72}
73
74template<typename T, unsigned D>
76{
77 return _cuboid.getDeltaR();
78}
79
80template<typename T, unsigned D>
82{
83 return _data[0][this->getCellId(latticeR)];
84}
85
86template<typename T, unsigned D>
87int BlockGeometry<T,D>::get(const int latticeR[D]) const
88{
89 return _data[0][this->getCellId(latticeR)];
90}
91
92template<typename T, unsigned D>
93int BlockGeometry<T,D>::get(std::size_t iCell) const
94{
95 return _data[0][iCell];
96}
98template<typename T, unsigned D>
101 if (this->isInside(latticeR)) {
102 return _data[0][this->getCellId(latticeR)];
103 }
104 else {
105 return 0;
106 }
107}
108
109template<typename T, unsigned D>
110int BlockGeometry<T,D>::getMaterial(const int latticeR[D]) const
112 LatticeR<D> loc{latticeR};
113 if (this->isInside(loc)) {
114 return _data[0][this->getCellId(loc)];
116 else {
117 return 0;
118 }
119}
120
121template<typename T, unsigned D>
122void BlockGeometry<T,D>::set(LatticeR<D> latticeR, int material)
123{
124 set(this->getCellId(latticeR), material);
127template<typename T, unsigned D>
128void BlockGeometry<T,D>::set(const int latticeR[D], int material)
129{
130 set(this->getCellId(latticeR), material);
131}
132
133template<typename T, unsigned D>
134void BlockGeometry<T,D>::set(std::size_t iCell, int material)
136 resetStatistics();
137 _data[0][iCell] = material;
138}
139
140template<typename T, unsigned D>
141void BlockGeometry<T,D>::getPhysR(T physR[D], const int latticeR[D]) const
142{
143 LatticeR<D> loc{latticeR};
144 getPhysR(physR, loc);
145 return;
146}
147
148template<typename T, unsigned D>
149void BlockGeometry<T,D>::getPhysR(T physR[D], LatticeR<D> loc) const
150{
151 _cuboid.getPhysR(physR, loc);
152 return;
153}
155template<typename T, unsigned D>
157{
158 return _iCglob;
159}
160
161template<typename T, unsigned D>
163{
164 if constexpr (D == 3) {
165 return Vector<int,3>(this->getNx(), this->getNy(), this->getNz());
166 } else {
167 return Vector<int,2>(this->getNx(), this->getNy());
168 }
169 __builtin_unreachable();
170}
171
172template<typename T, unsigned D>
173template<typename DESCRIPTOR>
174int BlockGeometry<T,D>::clean(bool verbose, std::vector<int> bulkMaterials)
175{
176 //using DESCRIPTOR = std::conditional_t<D==2,descriptors::D2Q5<>,descriptors::D3Q27<>>;
177int counter=0;
178 bool toClean = true;
179 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
180 // material not 0 and not in bulkMaterials
181 if (get(latticeR) != 0 && (! util::isContained(bulkMaterials, get(latticeR))) ) {
182 toClean = true;
183 for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
184 if ( util::isContained(bulkMaterials, get(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))) ){
185 toClean = false;
186 }
187 }
188 if (toClean){
189 set(latticeR, 0);
190 counter++;
192 }
193 });
194 if (verbose) {
195 clout << "cleaned "<< counter << " outer boundary voxel(s)" << std::endl;
196 }
197 return counter;
199
200template<typename T, unsigned D>
201int BlockGeometry<T,D>::outerClean(bool verbose, std::vector<int> bulkMaterials)
203 int counter=0;
204 using DESCRIPTOR = std::conditional_t<D==2,descriptors::D2Q9<>,descriptors::D3Q27<>>;
205 bool toClean = false;
206 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
207 if (util::isContained(bulkMaterials, get(latticeR))) {
208 toClean = false;
209 for(int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
210 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop))) == 0){
211 toClean = true;
212 }
213 }
214 if(toClean){
215 set(latticeR, 0);
216 counter++;
217 }
218 }
219 });
220 if (verbose) {
221 clout << "cleaned "<< counter << " outer fluid voxel(s)" << std::endl;
222 }
223 return counter;
224}
225
226template<typename T, unsigned D>
228{
229 int count2 = 0;
230 using DESCRIPTOR = std::conditional_t<D==2,descriptors::D2Q5<>,descriptors::D3Q7<>>;
231
232 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
233 if (get(latticeR) != 1 && get(latticeR) != 0) {
234 if constexpr (D==3){
235 bool var[7] = {false};
236 for(int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
237 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
238 var[iPop] = true;
239 }
240 }
241 int comb[12][3]={{1,4,2},{1,4,3},{1,4,5},{1,4,6},{2,5,1},{2,5,3},{2,5,4},{2,5,6},{3,6,1},{3,6,2},{3,6,4},{3,6,5}};
242 for(int i = 0; i < 12; i++){
243 if (var[(comb[i][0])] == true
244 && var[(comb[i][1])] == true
245 && var[(comb[i][2])] == true){
246 set(latticeR, 1);
247 count2++;//count2 is the same value as before, count2 gets increased even if this cell has already been cleaned
248 }
249 }
250 } else {
251 int var = 0;
252 for(int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
253 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
254 var = var+1;
255 }
256 }
257 if(var >= 3){
258 set(latticeR, 1);
259 count2++;//count2 differs from original count2
260 }
261 }
262 }
263 });
264
265 if (verbose) {
266 clout << "cleaned "<< count2 << " inner boundary voxel(s)" << std::endl;
267 }
268 return count2;
269}
270
271template<typename T, unsigned D>
272int BlockGeometry<T,D>::innerClean(int fromM, bool verbose)
273{
274 int count2 = 0;
275 using DESCRIPTOR = std::conditional_t<D==2,descriptors::D2Q5<>,descriptors::D3Q7<>>;
276
277 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
278 if (get(latticeR) != 1 && get(latticeR) != 0 && get(latticeR) == fromM) {
279 if constexpr (D==3){
280 bool var[7] = {false};
281 for(int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
282 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
283 var[iPop] = true;
284 }
285 }
286 int comb[12][3]={{1,4,2},{1,4,3},{1,4,5},{1,4,6},{2,5,1},{2,5,3},{2,5,4},{2,5,6},{3,6,1},{3,6,2},{3,6,4},{3,6,5}};
287 for(int i = 0; i < 12; i++){
288 if (var[(comb[i][0])] == true
289 && var[(comb[i][1])] == true
290 && var[(comb[i][2])] == true){
291 set(latticeR, 1);
292 count2++;//count2 is the same value as before, count2 gets increased even if this cell has already been cleaned
293 }
294 }
295 } else {
296 int var = 0;
297 for(int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
298 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
299 var = var+1;
300 }
301 }
302 if(var >= 3){
303 set(latticeR, 1);
304 count2++;//count2 differs from original count2
305 }
306 }
307 }
308 });
309
310 if (verbose){
311 clout << "cleaned "<< count2
312 << " inner boundary voxel(s) of Type " << fromM << std::endl;
313 }
314 return count2;
315}
316
317template<typename T, unsigned D>
319{
320 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
321 T physR[D] { };
322 bool output{};
323 getPhysR(physR, latticeR);
324 domain(&output, physR);
325 if (output) {
326 set(latticeR, 0);
327 }
328 });
329}
330
331template<typename T, unsigned D>
332bool BlockGeometry<T,D>::find(int material, std::vector<unsigned> offset,
333 std::vector<int> var)
334{
335 bool found = false;
336 for (var[0] = 0; var[0] < this->getNx(); var[0]++) {
337 for (var[1] = 0; var[1] < this->getNy(); var[1]++) {
338 if constexpr(D==3){
339 for (var[2] = 0; var[2] < this->getNz(); var[2]++) {
340 found = check(material, var, offset);
341 }
342 } else {
343 found = check(material, var, offset);
344 }
345 if (found) {
346 return found;
347 }
348 }
349 }
350 return found;
351}
352
353template<typename T, unsigned D>
354bool BlockGeometry<T,D>::check(int material, std::vector<int> var,
355 std::vector<unsigned> offset)
356{
357 bool found = true;
358 for (int iOffsetX = -offset[0]; iOffsetX <= (int) offset[0]; ++iOffsetX) {
359 for (int iOffsetY = -offset[1]; iOffsetY <= (int) offset[1]; ++iOffsetY) {
360 if constexpr (D==3){
361 for (int iOffsetZ = -offset[2]; iOffsetZ <= (int) offset[2]; ++iOffsetZ) {
362 if (getMaterial({var[0] + iOffsetX, var[1] + iOffsetY, var[2] + iOffsetZ}) != material) {
363 found = false;
364 }
365 }
366 } else {
367 if (getMaterial({var[0] + iOffsetX, var[1] + iOffsetY}) != material) {
368 found = false;
369 }
370 }
371 }
372 }
373 return found;
374}
375
376template<typename T, unsigned D>
378{
379 bool error = false;
380 using DESCRIPTOR = std::conditional_t<D==2,descriptors::D2Q9<>,descriptors::D3Q27<>>;
381 bool errorFound = false;
382 this->forSpatialLocations([&](LatticeR<D> latticeR)
383 {
384 if (get(latticeR) == 0) {
385 errorFound = false;
386 for(int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
387 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop))) == 1){
388 errorFound = true;
389 }
390 }
391 if(errorFound){
392 error = true;
393 }
394 }
395 });
396
397 if (verbose) {
398 if (error) {
399 clout << "error!" << std::endl;
400 }
401 else {
402 clout << "the model is correct!" << std::endl;
403 }
404 }
405 return error;
406}
407
408template<typename T, unsigned D>
409void BlockGeometry<T,D>::rename(int fromM, int toM)
410{
411 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
412 if (get(latticeR) == fromM) {
413 set(latticeR, toM);
414 }
415 });
416}
417
418template<typename T, unsigned D>
419void BlockGeometry<T,D>::rename(int fromM, int toM, IndicatorF<T,D>& condition)
420{
421 T physR[D];
422 this->forSpatialLocations([&](LatticeR<D> latticeR) {
423 if (get(latticeR) == fromM) {
424 getPhysR(physR, latticeR);
425 bool inside[1];
426 condition(inside, physR);
427 if (inside[0]) {
428 set(latticeR, toM);
429 }
430 }
431 });
432}
433
434template<typename T, unsigned D>
435void BlockGeometry<T,D>::rename(int fromM, int toM, LatticeR<D> offset)
436{
437 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
438 if (get(latticeR) == fromM) {
439 bool found = true;
440 for (int iOffsetX = -offset[0]; iOffsetX <= (int) offset[0]; ++iOffsetX) {
441 for (int iOffsetY = -offset[1]; iOffsetY <= (int) offset[1]; ++iOffsetY) {
442 if constexpr (D == 3) {
443 for (int iOffsetZ = -offset[2]; iOffsetZ <= (int) offset[2]; ++iOffsetZ) {
444 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY, latticeR[2] + iOffsetZ}) != fromM) {
445 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY, latticeR[2] + iOffsetZ}) != 1245) {
446 found = false;
447 }
448 }
449 }
450 } else {
451 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY}) != fromM) {
452 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY}) != 1245) {
453 found = false;
454 }
455 }
456 }
457 }
458 }
459 if (found) {
460 set(latticeR, 1245); // TODO Replace this embarassing approach
461 }
462 }
463 });
464 rename(1245,toM);
465}
466
467template<typename T, unsigned D>
468void BlockGeometry<T,D>::rename(int fromM, int toM, int testM,
469 std::vector<int> testDirection)
470{
471 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
472 if (get(latticeR) == fromM) {
473 // flag that indicates the renaming of the current voxel, valid voxels are not renamed
474 bool isValid = true;
475 for (int iOffsetX = util::min(testDirection[0],0); iOffsetX <= util::max(testDirection[0],0); ++iOffsetX) {
476 for (int iOffsetY = util::min(testDirection[1],0); iOffsetY <= util::max(testDirection[1],0); ++iOffsetY) {
477 if constexpr (D == 3){
478 for (int iOffsetZ = util::min(testDirection[2],0); iOffsetZ <= util::max(testDirection[2],0); ++iOffsetZ) {
479 if (iOffsetX!=0 || iOffsetY!=0 || iOffsetZ!=0) {
480 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY, latticeR[2] + iOffsetZ}) != testM) {
481 isValid = false;
482 }
483 }
484 }
485 } else {
486 if (iOffsetX!=0 || iOffsetY!=0) {
487 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY}) != testM) {
488 isValid = false;
489 }
490 }
491 }
492 }
493 }
494 if (!isValid) {
495 set(latticeR, toM);
496 }
497 }
498 });
499}
500
501
502template<typename T, unsigned D>
503void BlockGeometry<T,D>::rename(int fromM, int toM, int fluidM,
504 IndicatorF<T,D>& condition, Vector<int,D> discreteNormal)
505{
506 rename(fromM, toM, condition);
507 Vector<int,D> testDirection(discreteNormal);
508 T physR[D];
509 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
510 if (get(latticeR) == toM) {
511 getPhysR(physR, latticeR);
512 bool inside[1];
513 condition(inside, physR);
514 if (inside[0]) {
515 if constexpr (D==3){
516 if (getMaterial({latticeR[0]+testDirection[0],latticeR[1]+testDirection[1],latticeR[2]+testDirection[2]})!=fluidM ||
517 getMaterial({latticeR[0]+2*testDirection[0],latticeR[1]+2*testDirection[1],latticeR[2]+2*testDirection[2]})!=fluidM ||
518 getMaterial({latticeR[0]-testDirection[0],latticeR[1]-testDirection[1],latticeR[2]-testDirection[2]})!=0 ) {
519 set(latticeR, fromM);
520 }
521 } else {
522 if (getMaterial({latticeR[0]+testDirection[0],latticeR[1]+testDirection[1]}) != fluidM ||
523 getMaterial({latticeR[0]+2*testDirection[0],latticeR[1]+2*testDirection[1]}) != fluidM ||
524 getMaterial({latticeR[0]-testDirection[0],latticeR[1]-testDirection[1]}) != 0) {
525 set(latticeR, fromM);
526 }
527 }
528 }
529 }
530 });
531}
532
533template<typename T, unsigned D>
534void BlockGeometry<T,D>::rename(int fromM, int toM, int fluidM,
535 IndicatorF<T,D>& condition)
536{
537 rename(fromM, toM, condition);
538 std::vector<int> testDirection = getStatistics().computeDiscreteNormal(toM);
539 T physR[D];
540 // values that have been incorrectly changed from "fromM to "toM" get assigned back to "fromM"
541 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
542 if (get(latticeR) == toM) {
543 getPhysR(physR, latticeR);
544 bool inside[1];
545 condition(inside, physR);
546 if (inside[0]) {
547 if constexpr(D==3){
548 if (getMaterial({latticeR[0]+testDirection[0],latticeR[1]+testDirection[1],latticeR[2]+testDirection[2]})!=fluidM ||
549 getMaterial({latticeR[0]+2*testDirection[0],latticeR[1]+2*testDirection[1],latticeR[2]+2*testDirection[2]})!=fluidM ||
550 getMaterial({latticeR[0]-testDirection[0],latticeR[1]-testDirection[1],latticeR[2]-testDirection[2]})!=0 ) {
551 set(latticeR, fromM);
552 }
553 } else {
554 if (getMaterial({latticeR[0]+testDirection[0],latticeR[1]+testDirection[1]})!=fluidM ||
555 getMaterial({latticeR[0]+2*testDirection[0],latticeR[1]+2*testDirection[1]})!=fluidM ||
556 getMaterial({latticeR[0]-testDirection[0],latticeR[1]-testDirection[1]})!=0 ) {
557 set(latticeR, fromM);
558 }
559 }
560 }
561 }
562 });
563}
564
565template<typename T, unsigned D>
566void BlockGeometry<T,D>::copyMaterialLayer(IndicatorF3D<T>& condition, int discreteNormal[D], int numberOfLayers)
567{
568 T physR[D];
569 this->forCoreSpatialLocations([&](LatticeR<D> latticeR) {
570 getPhysR(physR, latticeR);
571 bool inside[1];
572 condition(inside, physR);
573 if (inside[0]) {
574 for (int i = 0; i < numberOfLayers; i++) {
575 if (0 <= latticeR[0] + i * discreteNormal[0] && latticeR[0] + i * discreteNormal[0] < this->getNx() &&
576 0 <= latticeR[1] + i * discreteNormal[1] && latticeR[1] + i * discreteNormal[1] < this->getNy()){
577 if constexpr(D==3){
578 if(0 <= latticeR[2] + i * discreteNormal[2] && latticeR[2] + i * discreteNormal[2] < this->getNz()){
579 set({latticeR[0] + i * discreteNormal[0], latticeR[1] + i * discreteNormal[1], latticeR[2] + i * discreteNormal[2]}, get(latticeR));
580 }
581 } else {
582 set({latticeR[0] + i * discreteNormal[0], latticeR[1] + i * discreteNormal[1]}, get(latticeR));
583 }
584 }
585 }
586 }
587 });
588}
589
590template<typename T, unsigned D>
591void BlockGeometry<T,D>::regionGrowing(int fromM, int toM, LatticeR<D> seed,
592 std::vector<int> offset, std::map<std::vector<int>, int>* tmp)
593{
594 std::map<std::vector<int>, int> tmp2;
595 bool firstCall = false;
596 if (tmp == nullptr) {
597 tmp = &tmp2;
598 firstCall = true;
599 }
600
601 if (getMaterial(seed) == fromM) {
602 std::vector<int> found;
603 found.push_back(seed[0]);
604 found.push_back(seed[1]);
605 if constexpr(D==3){
606 found.push_back(seed[2]);
607 if (tmp->count(found) == 0) {
608 (*tmp)[found] = 2;
609 if (offset[0] != 0) {
610 regionGrowing(fromM, toM, {seed[0] + 1, seed[1], seed[2]}, offset, tmp);
611 regionGrowing(fromM, toM, {seed[0] - 1, seed[1], seed[2]}, offset, tmp);
612 }
613 if (offset[1] != 0) {
614 regionGrowing(fromM, toM, {seed[0], seed[1] + 1, seed[2]}, offset, tmp);
615 regionGrowing(fromM, toM, {seed[0], seed[1] - 1, seed[2]}, offset, tmp);
616 }
617 if (offset[2] != 0) {
618 regionGrowing(fromM, toM, {seed[0], seed[1], seed[2] + 1}, offset, tmp);
619 regionGrowing(fromM, toM, {seed[0], seed[1], seed[2] - 1}, offset, tmp);
620 }
621 }
622 } else {
623 if (tmp->count(found) == 0) {
624 (*tmp)[found] = 2;
625 if (offset[0] != 0) {
626 regionGrowing(fromM, toM, {seed[0] + 1, seed[1]}, offset, tmp);
627 regionGrowing(fromM, toM, {seed[0] - 1, seed[1]}, offset, tmp);
628 }
629 if (offset[1] != 0) {
630 regionGrowing(fromM, toM, {seed[0], seed[1] + 1}, offset, tmp);
631 regionGrowing(fromM, toM, {seed[0], seed[1] - 1}, offset, tmp);
632 }
633 }
634 }
635 }
636 if (firstCall) {
637 std::map<std::vector<int>, int>::iterator iter;
638 for (iter = tmp->begin(); iter != tmp->end(); iter++) {
639 if constexpr(D==3){
640 set((iter->first)[0],(iter->first)[1],(iter->first)[2], toM);
641 } else {
642 set((iter->first)[0],(iter->first)[1], toM);
643 }
644 }
645 }
646 return;
647}
648
649template<typename T, unsigned D>
650void BlockGeometry<T,D>::printLayer(std::vector<int> min, std::vector<int> max, bool linenumber)//each call of this function must be carefully changed(previous:(xo,x1,y0,y1),now:({x0,yo},{x1,z1}))
651{
652 for (int x = min[0]; x <= max[0]; x++) {
653 if (linenumber) {
654 clout << x << ": ";
655 }
656 for (int y = min[1]; y <= max[1]; y++) {
657 if constexpr (D==3){
658 for (int z = min[2]; z <= max[2]; z++) {
659 clout << getMaterial({x, y, z}) << " ";
660 }
661 if (max[1] - min[1] != 0 && max[2] - min[2] != 0) {
662 clout << std::endl;
663 }
664 } else {
665 clout << getMaterial({x, y}) << " ";
666 }
667 }
668 if (max[0] - min[0] != 0) {
669 clout << std::endl;
670 }
671 }
672 clout << std::endl;
673}
674
675template<typename T, unsigned D>
676void BlockGeometry<T,D>::printLayer(int direction, int layer, bool linenumber)
677{
678 assert(direction >= 0 && direction <= 2);
679 if constexpr(D==3){
680 switch (direction) {
681 case 0:
682 printLayer({layer, 0, 0},{layer, this->getNy() - 1, this->getNz() - 1}, linenumber);
683 break;
684 case 1:
685 printLayer({0, layer, 0}, {this->getNx() - 1, layer, this->getNz() - 1}, linenumber);
686 break;
687 case 2:
688 printLayer({0, 0, layer}, {this->getNx() - 1, this->getNy() - 1, layer}, linenumber);
689 break;
690 }
691 } else {
692 switch (direction) {
693 case 0:
694 printLayer({layer, 0}, {layer, this->getNy() - 1}, linenumber);
695 break;
696 case 1:
697 printLayer({0, layer}, {this->getNx() - 1, layer}, linenumber);
698 break;
699 }
700 }
701}
702
703template<typename T, unsigned D>
704void BlockGeometry<T,D>::printNode(std::vector<int> loc)
705{
706 for (int x = loc[0] - 1; x <= loc[0] + 1; x++) {
707 clout << "x=" << x << std::endl;
708 for (int y = loc[1] - 1; y <= loc[1] + 1; y++) {
709 if constexpr (D==3){
710 for (int z = loc[2] - 1; z <= loc[2] + 1; z++) {
711 clout << getMaterial({x, y, z}) << " ";
712 }
713 clout << std::endl;
714 } else {
715 clout << getMaterial({x, y}) << " ";
716 }
717 }
718 clout << std::endl;
719 }
720 clout << std::endl;
721}
722
723template<typename T, unsigned D>
725{
726 _statistics.getStatisticsStatus() = true;
727}
728
729template<typename T, unsigned D>
731{
732 return _data.getNblock();
733}
734
735template<typename T, unsigned D>
737{
738 return _data.getSerializableSize();
739}
740
741template<typename T, unsigned D>
742bool* BlockGeometry<T,D>::getBlock(std::size_t iBlock, std::size_t& sizeBlock, bool loadingMode)
743{
744 std::size_t currentBlock = 0;
745 bool* dataPtr = nullptr;
746
747 this->registerSerializableOfConstSize(iBlock, sizeBlock, currentBlock, dataPtr, _data, loadingMode);
748
749 return dataPtr;
750}
751
752} // namespace olb
753
754#endif
Representation of the 2D block geometry view – header file.
Representation of a block geometry.
int innerClean(bool verbose=true)
Changes all cell materials which are not 0 or 1 to 1 if there is a non robust constiallation.
void rename(int fromM, int toM)
Replaces all material numbers (fromM) to another (toM)
void reset(IndicatorF< T, D > &domain)
Resets all cell materials inside of a domain to 0.
const BlockGeometryStatistics< T, D > & getStatistics() const
Read only access to the associated block statistic.
Vector< T, D > getPhysR(LatticeR< D > latticeR)
int getMaterial(LatticeR< D > latticeR) const
returns the (iX,iY) entry in the 2D scalar field
int clean(bool verbose=true, std::vector< int > bulkMaterials={1})
Changes all cell materials which are not in bulkMaterials to 0 if there is no neighbour from bulkMate...
BlockGeometry(Cuboid< T, D > &cuboid, int padding, int iCglob=-1)
bool find(int material, std::vector< unsigned > offset, std::vector< int > var)
Returns the coordinates (iX,iY) of a voxel with a given material number (material) if there exists an...
int outerClean(bool verbose=true, std::vector< int > bulkMaterials={1})
Changes all cell materials from bulkMaterials to 0 if there is a neighbour with material 0.
bool check(int material, std::vector< int > var, std::vector< unsigned > offset)
Returns true if at position (iX,iY) and in a neighbourhood of size (offsetX,offsetY) only voxels with...
int getIcGlob() const
Read only access to the global iC number which is given !=-1 if the block geometries are part of a su...
void regionGrowing(int fromM, int toM, LatticeR< D > seed, std::vector< int > offset, std::map< std::vector< int >, int > *tmp=nullptr)
Replaces all material numbers (fromM) to another (toM) using a seed point and max....
std::size_t getNblock() const override
Number of data blocks for the serializable interface.
T getDeltaR() const
Read only access to the voxel size given in SI units (meter)
bool checkForErrors(bool verbose=true) const
Checks for errors (searches for all outer voxels (=0) with an inner voxel (=1) as a direct neighbour)
void copyMaterialLayer(IndicatorF3D< T > &condition, int discreteNormal[D], int numberOfLayers)
Copy a layer of material numbers inside an indicator in a discrete normal direction.
bool * getBlock(std::size_t iBlock, std::size_t &sizeBlock, bool loadingMode) override
Return a pointer to the memory of the current block and its size for the serializable interface.
Vector< T, D > getOrigin() const
Read only access to the origin position given in SI units (meter)
void printNode(std::vector< int > loc)
Prints a chosen node and its neighbourhood.
Vector< int, D > getExtent() const
Returns the extend of the block in lattice units.
void printLayer(std::vector< int > min, std::vector< int > max, bool linenumber=false)
Prints a chosen part of the block geometry.
void set(LatticeR< D > latticeR, int material)
Write access to a material number.
std::size_t getSerializableSize() const override
Binary size for the serializer.
std::enable_if_t< sizeof...(L)==D, int > get(L... latticeR) const
Read-only access to a material number.
Base of a regular block.
IndicatorF3D is an application from .
Plain old scalar vector.
Definition vector.h:47
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
bool isContained(const C &c, U object)
Check, if object is contained in iteratable container c.
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
Top level namespace for all of OpenLB.
std::conditional_t< D==2, BlockGeometryStatistics2D< T >, BlockGeometryStatistics3D< T > > BlockGeometryStatistics
Definition aliases.h:178
std::conditional_t< D==2, IndicatorF2D< T >, IndicatorF3D< T > > IndicatorF
Definition aliases.h:258
std::conditional_t< D==2, Cuboid2D< T >, Cuboid3D< T > > Cuboid
Definition aliases.h:37