OpenLB 1.7
Loading...
Searching...
No Matches
blockGeometryStatistics3D.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 * 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
28#ifndef BLOCK_GEOMETRY_STATISTICS_3D_HH
29#define BLOCK_GEOMETRY_STATISTICS_3D_HH
30
31#include <iostream>
32#include <math.h>
33#include <fstream>
34#include <sstream>
35#include "utilities/omath.h"
36
37#include "geometry/blockGeometry.h"
38#include "geometry/blockGeometryStatistics3D.h"
39
40namespace olb {
41
42template<typename T>
43BlockGeometryStatistics3D<T>::BlockGeometryStatistics3D(BlockGeometry<T,3>* blockGeometry)
44 : _blockGeometry(blockGeometry),
45 clout(std::cout,"BlockGeometryStatistics3D")
46{
47 _statisticsUpdateNeeded = true;
48}
49
50template<typename T>
51bool& BlockGeometryStatistics3D<T>::getStatisticsStatus()
52{
53 return _statisticsUpdateNeeded;
54}
55
56template<typename T>
57bool const & BlockGeometryStatistics3D<T>::getStatisticsStatus() const
58{
59 return _statisticsUpdateNeeded;
60}
61
62
63template<typename T>
64void BlockGeometryStatistics3D<T>::update(bool verbose)
65{
66 const_this = const_cast<const BlockGeometryStatistics3D<T>*>(this);
67
68 if (getStatisticsStatus() ) {
69 _material2n.clear();
70 _blockGeometry->forCoreSpatialLocations([&](auto iX, auto iY, auto iZ) {
71 takeStatistics(iX,iY,iZ);
72 });
73
74 _nMaterials=int();
75 std::map<int, std::size_t>::iterator iter;
76 for (iter = _material2n.begin(); iter != _material2n.end(); iter++) {
77 _nMaterials++;
78 }
79
80 if (verbose) {
81 clout << "updated" << std::endl;
82 }
83 getStatisticsStatus()=false;
84 }
85}
86
87
88template<typename T>
89int BlockGeometryStatistics3D<T>::getNmaterials()
90{
91 update();
92 return const_this->getNmaterials();
93}
94
95template<typename T>
96int BlockGeometryStatistics3D<T>::getNmaterials() const
97{
98 return _nMaterials;
99}
100
101template<typename T>
102std::size_t BlockGeometryStatistics3D<T>::getNvoxel(int material)
103{
104 update();
105 return const_this->getNvoxel(material);
106}
107
108template<typename T>
109std::size_t BlockGeometryStatistics3D<T>::getNvoxel(int material) const
110{
111 try {
112 return _material2n.at(material);
113 }
114 catch (std::out_of_range& ex) {
115 return 0;
116 }
117}
118
119template<typename T>
120std::map<int, std::size_t> BlockGeometryStatistics3D<T>::getMaterial2n()
121{
122 update();
123 return const_this->getMaterial2n();
124}
125
126template<typename T>
127std::map<int, std::size_t> BlockGeometryStatistics3D<T>::getMaterial2n() const
128{
129 return _material2n;
130}
131
132template<typename T>
133std::size_t BlockGeometryStatistics3D<T>::getNvoxel()
134{
135 update();
136 return const_this->getNvoxel();
137}
138
139template<typename T>
140std::size_t BlockGeometryStatistics3D<T>::getNvoxel() const
141{
142 std::size_t total = 0;
143 for (const auto& material : _material2n) {
144 total += material.second;
145 }
146 return total;
147}
148template<typename T>
149std::vector<int> BlockGeometryStatistics3D<T>::getMinLatticeR(int material)
150{
151 update();
152 return const_this->getMinLatticeR(material);
153}
154
155template<typename T>
156std::vector<int> BlockGeometryStatistics3D<T>::getMinLatticeR(int material) const
157{
158 try {
159 return _material2min.at(material);
160 }
161 catch (std::out_of_range& ex) {
162 std::vector<int> null;
163 return null;
164 }
165}
166
167template<typename T>
168std::vector<int> BlockGeometryStatistics3D<T>::getMaxLatticeR(int material)
169{
170 update();
171 return const_this->getMaxLatticeR(material);
172}
173
174template<typename T>
175std::vector<int> BlockGeometryStatistics3D<T>::getMaxLatticeR(int material) const
176{
177 try {
178 return _material2max.at(material);
179 }
180 catch (std::out_of_range& ex) {
181 std::vector<int> null;
182 return null;
183 }
184}
185
186template<typename T>
187std::vector<T> BlockGeometryStatistics3D<T>::getMinPhysR(int material) const
188{
189 std::vector<T> tmp(3,T());
190 _blockGeometry->getPhysR(&(tmp[0]), &(getMinLatticeR(material)[0]));
191 return tmp;
192}
193
194template<typename T>
195std::vector<T> BlockGeometryStatistics3D<T>::getMaxPhysR(int material) const
196{
197 std::vector<T> tmp(3,T());
198 _blockGeometry->getPhysR(&(tmp[0]), &(getMaxLatticeR(material)[0]));
199 return tmp;
200}
201
202template<typename T>
203std::vector<T> BlockGeometryStatistics3D<T>::getLatticeExtend(int material)
204{
205 update();
206 return const_this->getLatticeExtend(material);
207}
208
209template<typename T>
210std::vector<T> BlockGeometryStatistics3D<T>::getLatticeExtend(int material) const
211{
212 try {
213 std::vector<T> extend;
214 for (int iDim = 0; iDim < 3; iDim++) {
215 extend.push_back(_material2max.at(material)[iDim] - _material2min.at(material)[iDim]);
216 }
217 return extend;
218 }
219 catch (std::out_of_range& ex) {
220 std::vector<T, std::allocator<T>> null;
221 return null;
222 }
223}
224
225template<typename T>
226std::vector<T> BlockGeometryStatistics3D<T>::getPhysExtend(int material)
227{
228 update();
229 return const_this->getPhysExtend(material);
230}
231
232template<typename T>
233std::vector<T> BlockGeometryStatistics3D<T>::getPhysExtend(int material) const
234{
235 std::vector<T> extend;
236 for (int iDim = 0; iDim < 3; iDim++) {
237 extend.push_back(getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim]);
238 }
239 return extend;
240}
241
242template<typename T>
243std::vector<T> BlockGeometryStatistics3D<T>::getPhysRadius(int material)
244{
245 update();
246 return const_this->getPhysRadius(material);
247}
248
249template<typename T>
250std::vector<T> BlockGeometryStatistics3D<T>::getPhysRadius(int material) const
251{
252 std::vector<T> radius;
253 for (int iDim=0; iDim<3; iDim++) {
254 radius.push_back((getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim])/2.);
255 }
256 return radius;
257}
258
259template<typename T>
260std::vector<T> BlockGeometryStatistics3D<T>::getCenterPhysR(int material)
261{
262 update();
263 return const_this->getCenterPhysR(material);
264}
265
266template<typename T>
267std::vector<T> BlockGeometryStatistics3D<T>::getCenterPhysR(int material) const
268{
269 std::vector<T> center;
270 for (int iDim=0; iDim<3; iDim++) {
271 center.push_back(getMinPhysR(material)[iDim] + getPhysRadius(material)[iDim]);
272 }
273 return center;
274}
275
276template<typename T>
277std::vector<int> BlockGeometryStatistics3D<T>::getType(int iX, int iY, int iZ, bool anyNormal)
278{
279 update();
280 return const_this->getType(iX, iY, iZ, anyNormal);
281}
282
283template<typename T>
284std::vector<int> BlockGeometryStatistics3D<T>::getType(const int* input, bool anyNormal) const
285{
286 return const_this->getType(input[0], input[1], input[2], anyNormal);
287}
288
289template<typename T>
290std::vector<int> BlockGeometryStatistics3D<T>::getType(int iX, int iY, int iZ, bool anyNormal) const
291{
292 std::vector<int> discreteNormal(4, 0);
293 std::vector<int> discreteNormal2(4, 0);
294 std::vector<int> nullVector(4, 0);
295
296 if (_blockGeometry->getMaterial({iX, iY, iZ}) != 1
297 && _blockGeometry->getMaterial({iX, iY, iZ}) != 0) {
298
299 //boundary0N and boundary 0P
300 if ( _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
301 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
302 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
303 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
304 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
305 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
306 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
307 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
308
309 if (_blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1) {
310 if (discreteNormal == nullVector) {
311 discreteNormal[0] = 0;
312 discreteNormal[1] = -1;
313 discreteNormal[2] = 0;
314 discreteNormal[3] = 0;
315 }
316 else {
317 discreteNormal2[0] = 0;
318 discreteNormal2[1] = -1;
319 discreteNormal2[2] = 0;
320 discreteNormal2[3] = 0;
321 }
322 }
323
324 if (_blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1) {
325 if (discreteNormal == nullVector) {
326 discreteNormal[0] = 0;
327 discreteNormal[1] = 1;
328 discreteNormal[2] = 0;
329 discreteNormal[3] = 0;
330 }
331 else {
332 discreteNormal2[0] = 0;
333 discreteNormal2[1] = 1;
334 discreteNormal2[2] = 0;
335 discreteNormal2[3] = 0;
336 }
337 }
338 }
339
340 // boundary1N and boundary1P
341 if ( _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
342 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
343 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
344 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
345 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
346 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
347 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
348 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0) {
349
350 if (_blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1) {
351 if (discreteNormal == nullVector) {
352 discreteNormal[0] = 0;
353 discreteNormal[1] = 0;
354 discreteNormal[2] = -1;
355 discreteNormal[3] = 0;
356 }
357 else {
358 discreteNormal2[0] = 0;
359 discreteNormal2[1] = 0;
360 discreteNormal2[2] = -1;
361 discreteNormal2[3] = 0;
362 }
363 }
364
365 if (_blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1) {
366 if (discreteNormal == nullVector) {
367 discreteNormal[0] = 0;
368 discreteNormal[1] = 0;
369 discreteNormal[2] = 1;
370 discreteNormal[3] = 0;
371 }
372 else {
373 discreteNormal2[0] = 0;
374 discreteNormal2[1] = 0;
375 discreteNormal2[2] = 1;
376 discreteNormal2[3] = 0;
377 }
378 }
379 }
380
381 // boundary2N and boundary2P
382 if (_blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
383 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
384 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
385 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
386 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
387 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
388 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
389 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
390
391 if (_blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1) {
392 if (discreteNormal == nullVector) {
393 discreteNormal[0] = 0;
394 discreteNormal[1] = 0;
395 discreteNormal[2] = 0;
396 discreteNormal[3] = -1;
397 }
398 else {
399 discreteNormal2[0] = 0;
400 discreteNormal2[1] = 0;
401 discreteNormal2[2] = 0;
402 discreteNormal2[3] = -1;
403 }
404 }
405
406 if (_blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1) {
407 if (discreteNormal == nullVector) {
408 discreteNormal[0] = 0;
409 discreteNormal[1] = 0;
410 discreteNormal[2] = 0;
411 discreteNormal[3] = 1;
412 }
413 else {
414 discreteNormal2[0] = 0;
415 discreteNormal2[1] = 0;
416 discreteNormal2[2] = 0;
417 discreteNormal2[3] = 1;
418 }
419 }
420 }
421
422 // externalCornerNNN and externalCornerNPN
423 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
424 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
425 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
426 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
427 && _blockGeometry->getMaterial({iX + 1, iY, iZ + 1}) != 1
428 && _blockGeometry->getMaterial({iX + 1, iY, iZ + 1}) != 0) {
429
430 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
431 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
432 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 1
433 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 0
434 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 1
435 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 0
436 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ + 1}) == 1) {
437
438 if (discreteNormal == nullVector) {
439 discreteNormal[0] = 1;
440 discreteNormal[1] = -1;
441 discreteNormal[2] = -1;
442 discreteNormal[3] = -1;
443 }
444 else {
445 discreteNormal2[0] = 1;
446 discreteNormal2[1] = -1;
447 discreteNormal2[2] = -1;
448 discreteNormal2[3] = -1;
449 }
450 }
451
452 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
453 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
454 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 1
455 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 0
456 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 1
457 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 0
458 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ + 1}) == 1) {
459
460 if (discreteNormal == nullVector) {
461 discreteNormal[0] = 1;
462 discreteNormal[1] = -1;
463 discreteNormal[2] = 1;
464 discreteNormal[3] = -1;
465 }
466 else {
467 discreteNormal2[0] = 1;
468 discreteNormal2[1] = -1;
469 discreteNormal2[2] = 1;
470 discreteNormal2[3] = -1;
471 }
472 }
473 }
474
475 // externalCornerNPP and externalCornerNNP
476 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
477 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
478 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
479 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
480 && _blockGeometry->getMaterial({iX + 1, iY, iZ - 1}) != 1
481 && _blockGeometry->getMaterial({iX + 1, iY, iZ - 1}) != 0) {
482
483 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
484 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
485 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 1
486 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 0
487 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 1
488 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 0
489 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ - 1}) == 1) {
490
491 if (discreteNormal == nullVector) {
492 discreteNormal[0] = 1;
493 discreteNormal[1] = -1;
494 discreteNormal[2] = 1;
495 discreteNormal[3] = 1;
496 }
497 else {
498 discreteNormal2[0] = 1;
499 discreteNormal2[1] = -1;
500 discreteNormal2[2] = 1;
501 discreteNormal2[3] = 1;
502 }
503 }
504
505 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
506 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
507 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 1
508 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 0
509 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 1
510 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 0
511 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ - 1}) == 1) {
512
513 if (discreteNormal == nullVector) {
514 discreteNormal[0] = 1;
515 discreteNormal[1] = -1;
516 discreteNormal[2] = -1;
517 discreteNormal[3] = 1;
518 }
519
520 else {
521 discreteNormal2[0] = 1;
522 discreteNormal2[1] = -1;
523 discreteNormal2[2] = -1;
524 discreteNormal2[3] = 1;
525 }
526 }
527 }
528
529 // externalCornerPPP and externalCornerPNP
530 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
531 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
532 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
533 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
534 && _blockGeometry->getMaterial({iX - 1, iY, iZ - 1}) != 1
535 && _blockGeometry->getMaterial({iX - 1, iY, iZ - 1}) != 0) {
536
537 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
538 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
539 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 1
540 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 0
541 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 1
542 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 0
543 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ - 1}) == 1) {
544
545 if (discreteNormal == nullVector) {
546 discreteNormal[0] = 1;
547 discreteNormal[1] = 1;
548 discreteNormal[2] = 1;
549 discreteNormal[3] = 1;
550 }
551 else {
552 discreteNormal2[0] = 1;
553 discreteNormal2[1] = 1;
554 discreteNormal2[2] = 1;
555 discreteNormal2[3] = 1;
556 }
557 }
558
559 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
560 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
561 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 1
562 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 0
563 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 1
564 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 0
565 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ - 1}) == 1) {
566
567 if (discreteNormal == nullVector) {
568 discreteNormal[0] = 1;
569 discreteNormal[1] = 1;
570 discreteNormal[2] = -1;
571 discreteNormal[3] = 1;
572 }
573 else {
574 discreteNormal2[0] = 1;
575 discreteNormal2[1] = 1;
576 discreteNormal2[2] = -1;
577 discreteNormal2[3] = 1;
578 }
579 }
580 }
581
582 // externalCornerPNN and externalCornerPPN
583 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
584 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
585 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
586 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
587 && _blockGeometry->getMaterial({iX - 1, iY, iZ + 1}) != 1
588 && _blockGeometry->getMaterial({iX - 1, iY, iZ + 1}) != 0) {
589
590 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
591 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
592 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 1
593 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 0
594 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 1
595 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 0
596 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ + 1}) == 1) {
597
598 if (discreteNormal == nullVector) {
599 discreteNormal[0] = 1;
600 discreteNormal[1] = 1;
601 discreteNormal[2] = -1;
602 discreteNormal[3] = -1;
603 }
604 else {
605 discreteNormal2[0] = 1;
606 discreteNormal2[1] = 1;
607 discreteNormal2[2] = -1;
608 discreteNormal2[3] = -1;
609 }
610 }
611
612 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
613 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
614 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 1
615 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 0
616 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 1
617 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 0
618 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ + 1}) == 1) {
619
620 if (discreteNormal == nullVector) {
621
622 discreteNormal[0] = 1;
623 discreteNormal[1] = 1;
624 discreteNormal[2] = 1;
625 discreteNormal[3] = -1;
626 }
627
628 else {
629 discreteNormal2[0] = 1;
630 discreteNormal2[1] = 1;
631 discreteNormal2[2] = 1;
632 discreteNormal2[3] = -1;
633 }
634 }
635 }
636
637 // internalCornerPPP and internalCornerPNP
638 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1
639 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1
640 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
641 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
642 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
643 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0) {
644
645 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
646 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
647 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
648
649 if (discreteNormal == nullVector) {
650 discreteNormal[0] = 2;
651 discreteNormal[1] = 1;
652 discreteNormal[2] = 1;
653 discreteNormal[3] = 1;
654 }
655 else {
656 discreteNormal2[0] = 2;
657 discreteNormal2[1] = 1;
658 discreteNormal2[2] = 1;
659 discreteNormal2[3] = 1;
660 }
661 }
662
663 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
664 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
665 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
666
667 if (discreteNormal == nullVector) {
668 discreteNormal[0] = 2;
669 discreteNormal[1] = 1;
670 discreteNormal[2] = -1;
671 discreteNormal[3] = 1;
672 }
673
674 else {
675 discreteNormal2[0] = 2;
676 discreteNormal2[1] = 1;
677 discreteNormal2[2] = -1;
678 discreteNormal2[3] = 1;
679 }
680 }
681 }
682
683 // internalCornerPNN and InternalCornerPPN
684 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1
685 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1
686 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
687 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
688 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
689 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0) {
690
691 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
692 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
693 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
694
695 if (discreteNormal == nullVector) {
696 discreteNormal[0] = 2;
697 discreteNormal[1] = 1;
698 discreteNormal[2] = -1;
699 discreteNormal[3] = -1;
700 }
701 else {
702 discreteNormal2[0] = 2;
703 discreteNormal2[1] = 1;
704 discreteNormal2[2] = -1;
705 discreteNormal2[3] = -1;
706 }
707 }
708
709 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
710 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
711 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
712
713 if (discreteNormal == nullVector) {
714 discreteNormal[0] = 2;
715 discreteNormal[1] = 1;
716 discreteNormal[2] = 1;
717 discreteNormal[3] = -1;
718 }
719 else {
720 discreteNormal2[0] = 2;
721 discreteNormal2[1] = 1;
722 discreteNormal2[2] = 1;
723 discreteNormal2[3] = -1;
724 }
725 }
726 }
727
728 // internalCornerNPP and internalCornerNNP
729 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1
730 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1
731 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
732 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
733 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
734 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0) {
735
736 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
737 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
738 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
739
740 if (discreteNormal == nullVector) {
741 discreteNormal[0] = 2;
742 discreteNormal[1] = -1;
743 discreteNormal[2] = 1;
744 discreteNormal[3] = 1;
745 }
746 else {
747 discreteNormal2[0] = 2;
748 discreteNormal2[1] = -1;
749 discreteNormal2[2] = 1;
750 discreteNormal2[3] = 1;
751 }
752 }
753
754 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
755 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
756 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
757
758 if (discreteNormal == nullVector) {
759 discreteNormal[0] = 2;
760 discreteNormal[1] = -1;
761 discreteNormal[2] = -1;
762 discreteNormal[3] = 1;
763 }
764
765 else {
766 discreteNormal2[0] = 2;
767 discreteNormal2[1] = -1;
768 discreteNormal2[2] = -1;
769 discreteNormal2[3] = 1;
770 }
771 }
772 }
773
774 // internalCornerNPN and internalCornerNNN
775 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1
776 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1
777 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
778 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
779 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
780 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0) {
781
782 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
783 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
784 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
785
786 if (discreteNormal == nullVector) {
787 discreteNormal[0] = 2;
788 discreteNormal[1] = -1;
789 discreteNormal[2] = 1;
790 discreteNormal[3] = -1;
791 }
792
793 else {
794 discreteNormal2[0] = 2;
795 discreteNormal2[1] = -1;
796 discreteNormal2[2] = 1;
797 discreteNormal2[3] = -1;
798 }
799 }
800
801 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
802 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
803 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
804
805 if (discreteNormal == nullVector) {
806 discreteNormal[0] = 2;
807 discreteNormal[1] = -1;
808 discreteNormal[2] = -1;
809 discreteNormal[3] = -1;
810 }
811 else {
812 discreteNormal2[0] = 2;
813 discreteNormal2[1] = -1;
814 discreteNormal2[2] = -1;
815 discreteNormal2[3] = -1;
816 }
817 }
818 }
819
820 // externalEdge0PN and externalEdge0NN
821 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
822 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
823 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
824 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
825 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
826 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
827 && _blockGeometry->getMaterial({iX + 1, iY, iZ + 1}) != 1
828 && _blockGeometry->getMaterial({iX - 1, iY, iZ + 1}) != 1) {
829
830 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) == 1
831 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
832 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
833
834 if (discreteNormal == nullVector) {
835 discreteNormal[0] = 3;
836 discreteNormal[1] = 0;
837 discreteNormal[2] = 1;
838 discreteNormal[3] = -1;
839 }
840 else {
841 discreteNormal2[0] = 3;
842 discreteNormal2[1] = 0;
843 discreteNormal2[2] = 1;
844 discreteNormal2[3] = -1;
845 }
846 }
847
848 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) == 1
849 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
850 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
851
852 if (discreteNormal == nullVector) {
853 discreteNormal[0] = 3;
854 discreteNormal[1] = 0;
855 discreteNormal[2] = -1;
856 discreteNormal[3] = -1;
857 }
858 else {
859 discreteNormal2[0] = 3;
860 discreteNormal2[1] = 0;
861 discreteNormal2[2] = -1;
862 discreteNormal2[3] = -1;
863 }
864 }
865 }
866
867 // externalEdge0NP and externalEdge0PP
868 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
869 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
870 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
871 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
872 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
873 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
874 && _blockGeometry->getMaterial({iX + 1, iY, iZ - 1}) != 1
875 && _blockGeometry->getMaterial({iX - 1, iY, iZ - 1}) != 1) {
876
877 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) == 1
878 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1) {
879 if (discreteNormal == nullVector) {
880 discreteNormal[0] = 3;
881 discreteNormal[1] = 0;
882 discreteNormal[2] = -1;
883 discreteNormal[3] = 1;
884 }
885 else {
886 discreteNormal2[0] = 3;
887 discreteNormal2[1] = 0;
888 discreteNormal2[2] = -1;
889 discreteNormal2[3] = 1;
890 }
891 }
892
893 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) == 1
894 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1) {
895 if (discreteNormal == nullVector) {
896 discreteNormal[0] = 3;
897 discreteNormal[1] = 0;
898 discreteNormal[2] = 1;
899 discreteNormal[3] = 1;
900 }
901 else {
902 discreteNormal2[0] = 3;
903 discreteNormal2[1] = 0;
904 discreteNormal2[2] = 1;
905 discreteNormal2[3] = 1;
906 }
907 }
908 }
909
910 // externalEdge1NN and externalEdge1NP
911 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
912 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
913 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
914 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
915 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
916 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0) {
917
918 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ + 1}) == 1
919 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
920 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
921 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
922 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1) {
923
924 if (discreteNormal == nullVector) {
925 discreteNormal[0] = 3;
926 discreteNormal[1] = -1;
927 discreteNormal[2] = 0;
928 discreteNormal[3] = -1;
929 }
930 else {
931 discreteNormal2[0] = 3;
932 discreteNormal2[1] = -1;
933 discreteNormal2[2] = 0;
934 discreteNormal2[3] = -1;
935 }
936 }
937
938 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ + 1}) == 1
939 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
940 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
941 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
942 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1) {
943
944 if (discreteNormal == nullVector) {
945 discreteNormal[0] = 3;
946 discreteNormal[1] = 1;
947 discreteNormal[2] = 0;
948 discreteNormal[3] = -1;
949 }
950 else {
951 discreteNormal2[0] = 3;
952 discreteNormal2[1] = 1;
953 discreteNormal2[2] = 0;
954 discreteNormal2[3] = -1;
955 }
956 }
957 }
958
959 // externalEdge1PN and externalEdge1PP
960 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
961 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
962 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
963 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
964 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
965 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0) {
966
967 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ - 1}) == 1
968 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
969 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
970 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
971 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1) {
972
973 if (discreteNormal == nullVector) {
974 discreteNormal[0] = 3;
975 discreteNormal[1] = -1;
976 discreteNormal[2] = 0;
977 discreteNormal[3] = 1;
978 }
979 else {
980 discreteNormal2[0] = 3;
981 discreteNormal2[1] = -1;
982 discreteNormal2[2] = 0;
983 discreteNormal2[3] = 1;
984 }
985 }
986
987 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ - 1}) == 1
988 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
989 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
990 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
991 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1) {
992
993 if (discreteNormal == nullVector) {
994 discreteNormal[0] = 3;
995 discreteNormal[1] = 1;
996 discreteNormal[2] = 0;
997 discreteNormal[3] = 1;
998 }
999 else {
1000 discreteNormal2[0] = 3;
1001 discreteNormal2[1] = 1;
1002 discreteNormal2[2] = 0;
1003 discreteNormal2[3] = 1;
1004 }
1005 }
1006 }
1007
1008 // externalEdge2NN and externalEdge2PN
1009 if ( _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1010 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1011 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1012 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1013 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
1014 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
1015
1016 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
1017 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) == 1) {
1018 if (discreteNormal == nullVector) {
1019 discreteNormal[0] = 3;
1020 discreteNormal[1] = -1;
1021 discreteNormal[2] = -1;
1022 discreteNormal[3] = 0;
1023 }
1024 else {
1025 discreteNormal2[0] = 3;
1026 discreteNormal2[1] = -1;
1027 discreteNormal2[2] = -1;
1028 discreteNormal2[3] = 0;
1029 }
1030 }
1031
1032 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
1033 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) == 1) {
1034 if (discreteNormal == nullVector) {
1035 discreteNormal[0] = 3;
1036 discreteNormal[1] = 1;
1037 discreteNormal[2] = -1;
1038 discreteNormal[3] = 0;
1039 }
1040 else {
1041 discreteNormal2[0] = 3;
1042 discreteNormal2[1] = 1;
1043 discreteNormal2[2] = -1;
1044 discreteNormal2[3] = 0;
1045 }
1046 }
1047 }
1048
1049 // externalEdge2PP and externalEdge2NP
1050 if ( _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1051 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1052 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1053 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1054 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
1055 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
1056
1057 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
1058 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) == 1) {
1059 if (discreteNormal == nullVector) {
1060 discreteNormal[0] = 3;
1061 discreteNormal[1] = 1;
1062 discreteNormal[2] = 1;
1063 discreteNormal[3] = 0;
1064 }
1065 else {
1066 discreteNormal2[0] = 3;
1067 discreteNormal2[1] = 1;
1068 discreteNormal2[2] = 1;
1069 discreteNormal2[3] = 0;
1070 }
1071 }
1072
1073 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
1074 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) == 1) {
1075 if (discreteNormal == nullVector) {
1076 discreteNormal[0] = 3;
1077 discreteNormal[1] = -1;
1078 discreteNormal[2] = 1;
1079 discreteNormal[3] = 0;
1080 }
1081 else {
1082 discreteNormal2[0] = 3;
1083 discreteNormal2[1] = -1;
1084 discreteNormal2[2] = 1;
1085 discreteNormal2[3] = 0;
1086 }
1087 }
1088 }
1089
1090 // internalEdge0NN and internalEdge0PN
1091 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
1092 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
1093 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
1094 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
1095 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1096 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1097 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1) {
1098
1099 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
1100 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
1101 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
1102 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 1
1103 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 0
1104 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 1
1105 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 0) {
1106
1107 discreteNormal[0] = 4;
1108 discreteNormal[1] = 0;
1109 discreteNormal[2] = -1;
1110 discreteNormal[3] = -1;
1111 }
1112 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
1113 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
1114 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
1115 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 1
1116 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 0
1117 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 1
1118 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 0) {
1119
1120 discreteNormal[0] = 4;
1121 discreteNormal[1] = 0;
1122 discreteNormal[2] = 1;
1123 discreteNormal[3] = -1;
1124 }
1125 }
1126
1127 // internalEdge0NP and internalEdge0PP
1128 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
1129 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
1130 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
1131 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
1132 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1133 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1134 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1) {
1135
1136 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
1137 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 1
1138 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) != 0
1139 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
1140 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
1141 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 1
1142 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) != 0) {
1143
1144 discreteNormal[0] = 4;
1145 discreteNormal[1] = 0;
1146 discreteNormal[2] = -1;
1147 discreteNormal[3] = 1;
1148 }
1149
1150 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
1151 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 1
1152 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) != 0
1153 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
1154 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
1155 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 1
1156 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) != 0) {
1157
1158 discreteNormal[0] = 4;
1159 discreteNormal[1] = 0;
1160 discreteNormal[2] = 1;
1161 discreteNormal[3] = 1;
1162 }
1163 }
1164
1165 // internalEdge1PP and internalEdge 1NP
1166 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1
1167 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 1
1168 && _blockGeometry->getMaterial({iX + 1, iY, iZ}) != 0
1169 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
1170 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0
1171 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
1172 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
1173
1174 if ( _blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1
1175 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1176 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1177 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 1
1178 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 0
1179 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 1
1180 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 0) {
1181
1182 discreteNormal[0] = 4;
1183 discreteNormal[1] = 1;
1184 discreteNormal[2] = 0;
1185 discreteNormal[3] = 1;
1186 }
1187
1188 if ( _blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1
1189 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1190 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1191 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 1
1192 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 0
1193 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 1
1194 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 0) {
1195
1196 discreteNormal[0] = 4;
1197 discreteNormal[1] = 1;
1198 discreteNormal[2] = 0;
1199 discreteNormal[3] = -1;
1200 }
1201 }
1202
1203 // internalEdge1PN and internalEdge1NN
1204 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1
1205 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 1
1206 && _blockGeometry->getMaterial({iX - 1, iY, iZ}) != 0
1207 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
1208 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0
1209 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
1210 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
1211
1212 if ( _blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1
1213 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1214 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1215 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 1
1216 && _blockGeometry->getMaterial({iX, iY + 1, iZ + 1}) != 0
1217 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 1
1218 && _blockGeometry->getMaterial({iX, iY - 1, iZ + 1}) != 0) {
1219
1220 discreteNormal[0] = 4;
1221 discreteNormal[1] = -1;
1222 discreteNormal[2] = 0;
1223 discreteNormal[3] = 1;
1224 }
1225
1226 if ( _blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1
1227 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1228 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1229 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 1
1230 && _blockGeometry->getMaterial({iX, iY + 1, iZ - 1}) != 0
1231 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 1
1232 && _blockGeometry->getMaterial({iX, iY - 1, iZ - 1}) != 0) {
1233
1234 discreteNormal[0] = 4;
1235 discreteNormal[1] = -1;
1236 discreteNormal[2] = 0;
1237 discreteNormal[3] = -1;
1238 }
1239 }
1240
1241 // internalEdge2PP and internalEdge2NP
1242 if ( _blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1
1243 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1244 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1245 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1246 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1247 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 1
1248 && _blockGeometry->getMaterial({iX, iY + 1, iZ}) != 0) {
1249
1250 if (_blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1
1251 && _blockGeometry->getMaterial({iX - 1, iY - 1, iZ}) == 1) {
1252
1253 discreteNormal[0] = 4;
1254 discreteNormal[1] = 1;
1255 discreteNormal[2] = 1;
1256 discreteNormal[3] = 0;
1257 }
1258
1259 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1
1260 && _blockGeometry->getMaterial({iX + 1, iY - 1, iZ}) == 1) {
1261
1262 discreteNormal[0] = 4;
1263 discreteNormal[1] = -1;
1264 discreteNormal[2] = 1;
1265 discreteNormal[3] = 0;
1266 }
1267 }
1268
1269 // internalEdge2PN and internalEdge2NN
1270 if ( _blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1
1271 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 1
1272 && _blockGeometry->getMaterial({iX, iY, iZ - 1}) != 0
1273 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 1
1274 && _blockGeometry->getMaterial({iX, iY, iZ + 1}) != 0
1275 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 1
1276 && _blockGeometry->getMaterial({iX, iY - 1, iZ}) != 0) {
1277
1278 if ( _blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1
1279 && _blockGeometry->getMaterial({iX - 1, iY + 1, iZ}) == 1) {
1280
1281 discreteNormal[0] = 4;
1282 discreteNormal[1] = 1;
1283 discreteNormal[2] = -1;
1284 discreteNormal[3] = 0;
1285 }
1286
1287 if ( _blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1
1288 && _blockGeometry->getMaterial({iX + 1, iY + 1, iZ}) == 1) {
1289
1290 discreteNormal[0] = 4;
1291 discreteNormal[1] = -1;
1292 discreteNormal[2] = -1;
1293 discreteNormal[3] = 0;
1294 }
1295 }
1296
1297 // special boundary conditions
1298 if (discreteNormal2 != nullVector && anyNormal == false) {
1299 discreteNormal = checkExtraBoundary(discreteNormal, discreteNormal2);
1300 }
1301 }
1302
1303 if (discreteNormal[1] == 0 && discreteNormal[2] == 0 && discreteNormal[3] == 0) {
1304 clout << "WARNING: no discreteNormal is found" << std::endl;
1305 }
1306 else if (_blockGeometry->getMaterial({iX-discreteNormal[1], iY-discreteNormal[2], iZ-discreteNormal[3]}) != 1) {
1307#ifdef OLB_DEBUG
1308 clout << "WARNING: discreteNormal is not pointing outside the fluid. Use option: anyNormal" << std::endl;
1309#endif
1310 }
1311
1312 return discreteNormal;
1313}
1314
1315template<typename T>
1316std::vector<int> BlockGeometryStatistics3D<T>::computeNormal(int iX, int iY, int iZ)
1317{
1318 update();
1319 return const_this->computeNormal(iX, iY, iZ);
1320}
1321
1322template<typename T>
1323std::vector<int> BlockGeometryStatistics3D<T>::computeNormal(int iX, int iY, int iZ) const
1324{
1325 std::vector<int> normal (3,int(0));
1326
1327 if (iX != 0) {
1328 if (_blockGeometry->getMaterial({iX - 1, iY, iZ}) == 1) {
1329 normal[0] = -1;
1330 }
1331 }
1332 if (iX != _nX - 1) {
1333 if (_blockGeometry->getMaterial({iX + 1, iY, iZ}) == 1) {
1334 normal[0] = 1;
1335 }
1336 }
1337 if (iY != 0) {
1338 if (_blockGeometry->getMaterial({iX, iY - 1, iZ}) == 1) {
1339 normal[1] = -1;
1340 }
1341 }
1342 if (iY != _nY - 1) {
1343 if (_blockGeometry->getMaterial({iX, iY + 1, iZ}) == 1) {
1344 normal[1] = 1;
1345 }
1346 }
1347 if (iZ != 0) {
1348 if (_blockGeometry->getMaterial({iX, iY, iZ - 1}) == 1) {
1349 normal[2] = -1;
1350 }
1351 }
1352 if (iZ != _nZ - 1) {
1353 if (_blockGeometry->getMaterial({iX, iY, iZ + 1}) == 1) {
1354 normal[2] = 1;
1355 }
1356 }
1357 return normal;
1358}
1359
1360template<typename T>
1361std::vector<T> BlockGeometryStatistics3D<T>::computeNormal(int material)
1362{
1363
1364 update();
1365 return const_this->computeNormal(material);
1366}
1367
1368template<typename T>
1369std::vector<T> BlockGeometryStatistics3D<T>::computeNormal(int material) const
1370{
1371 std::vector<T> normal (3,int(0));
1372 std::vector<int> minC = getMinLatticeR(material);
1373 std::vector<int> maxC = getMaxLatticeR(material);
1374 for (int iX = minC[0]; iX<=maxC[0]; iX++) {
1375 for (int iY = minC[1]; iY<=maxC[1]; iY++) {
1376 for (int iZ = minC[2]; iZ<=maxC[2]; iZ++) {
1377 if (_blockGeometry->getMaterial({iX,iY,iZ}) == material) {
1378 auto n = computeNormal(iX,iY,iZ);
1379 normal[0]+=n[0];
1380 normal[1]+=n[1];
1381 normal[2]+=n[2];
1382 }
1383 }
1384 }
1385 }
1386 T norm = util::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
1387 if (norm>0.) {
1388 normal[0]/=norm;
1389 normal[1]/=norm;
1390 normal[2]/=norm;
1391 }
1392 return normal;
1393}
1394
1395template<typename T>
1396std::vector<int> BlockGeometryStatistics3D<T>::computeDiscreteNormal(int material, T maxNorm)
1397{
1398 update();
1399 return const_this->computeDiscreteNormal(material, maxNorm);
1400}
1401
1402template<typename T>
1403std::vector<int> BlockGeometryStatistics3D<T>::computeDiscreteNormal(int material, T maxNorm) const
1404{
1405 std::vector<T> normal = computeNormal(material);
1406 std::vector<int> discreteNormal(3,int(0));
1407
1408 T smallestAngle = T(0);
1409 for (int iX = -1; iX<=1; iX++) {
1410 for (int iY = -1; iY<=1; iY++) {
1411 for (int iZ = -1; iZ<=1; iZ++) {
1412 T norm = util::sqrt(iX*iX+iY*iY+iZ*iZ);
1413 if (norm>0.&& norm<maxNorm) {
1414 T angle = (iX*normal[0] + iY*normal[1] + iZ*normal[2])/norm;
1415 if (angle>=smallestAngle) {
1416 smallestAngle=angle;
1417 discreteNormal[0] = iX;
1418 discreteNormal[1] = iY;
1419 discreteNormal[2] = iZ;
1420 }
1421 }
1422 }
1423 }
1424 }
1425 return discreteNormal;
1426}
1427
1428template<typename T>
1429bool BlockGeometryStatistics3D<T>::check(int material, int iX, int iY,
1430 int iZ, unsigned offsetX, unsigned offsetY, unsigned offsetZ)
1431{
1432 update();
1433 return const_this->check(material, iX, iY, iZ, offsetX, offsetY, offsetZ);
1434}
1435
1436template<typename T>
1437bool BlockGeometryStatistics3D<T>::check(int material, int iX, int iY,
1438 int iZ, unsigned offsetX, unsigned offsetY, unsigned offsetZ) const
1439{
1440 bool found = true;
1441 for (int iOffsetX = -offsetX; iOffsetX <= (int) offsetX; ++iOffsetX) {
1442 for (int iOffsetY = -offsetY; iOffsetY <= (int) offsetY; ++iOffsetY) {
1443 for (int iOffsetZ = -offsetZ; iOffsetZ <= (int) offsetZ; ++iOffsetZ) {
1444 if (_blockGeometry->getMaterial({iX + iOffsetX, iY + iOffsetY,
1445 iZ + iOffsetZ}) != material) {
1446 found = false;
1447 }
1448 }
1449 }
1450 }
1451 return found;
1452}
1453
1454template<typename T>
1455bool BlockGeometryStatistics3D<T>::find(int material, unsigned offsetX,
1456 unsigned offsetY, unsigned offsetZ, int& foundX, int& foundY,
1457 int& foundZ)
1458{
1459 update();
1460 return const_this->find(material, offsetX, offsetY, offsetZ, foundX, foundY, foundZ);
1461}
1462
1463template<typename T>
1464bool BlockGeometryStatistics3D<T>::find(int material, unsigned offsetX,
1465 unsigned offsetY, unsigned offsetZ, int& foundX, int& foundY,
1466 int& foundZ) const
1467{
1468 bool found = false;
1469 for (foundX = 0; foundX < _nX; foundX++) {
1470 for (foundY = 0; foundY < _nY; foundY++) {
1471 for (foundZ = 0; foundZ < _nZ; foundZ++) {
1472 found = check(material, foundX, foundY, foundZ, offsetX,
1473 offsetY, offsetZ);
1474 if (found) {
1475 return found;
1476 }
1477 }
1478 }
1479 }
1480 return found;
1481}
1482
1483template<typename T>
1484void BlockGeometryStatistics3D<T>::print()
1485{
1486 update();
1487 return const_this->print();
1488}
1489
1490template<typename T>
1491void BlockGeometryStatistics3D<T>::print() const
1492{
1493 try {
1494 for (const auto& material : _material2n) {
1495 clout << "materialNumber=" << material.first
1496 << "; count=" << material.second
1497 << "; minLatticeR=(" << _material2min.at(material.first)[0] <<","
1498 << _material2min.at(material.first)[1] <<","<< _material2min.at(material.first)[2] <<")"
1499 << "; maxLatticeR=(" << _material2max.at(material.first)[0] <<","
1500 << _material2max.at(material.first)[1] <<","<< _material2max.at(material.first)[2] <<")"
1501 << std::endl;
1502 }
1503 }
1504 catch (std::out_of_range& ex) {
1505 }
1506}
1507
1508template<typename T>
1509void BlockGeometryStatistics3D<T>::takeStatistics(int iX, int iY, int iZ)
1510{
1511 int type = _blockGeometry->getMaterial({iX, iY, iZ});
1512 if (_material2n.count(type) == 0) {
1513 _material2n[type] = 1;
1514 std::vector<int> minCo;
1515 std::vector<int> maxCo;
1516 minCo.push_back(iX);
1517 minCo.push_back(iY);
1518 minCo.push_back(iZ);
1519 _material2min[type] = minCo;
1520 maxCo.push_back(iX);
1521 maxCo.push_back(iY);
1522 maxCo.push_back(iZ);
1523 _material2max[type] = maxCo;
1524
1525 }
1526 else {
1527 _material2n[type]++;
1528 if (iX < _material2min[type][0]) {
1529 _material2min[type][0] = iX;
1530 }
1531 if (iY < _material2min[type][1]) {
1532 _material2min[type][1] = iY;
1533 }
1534 if (iZ < _material2min[type][2]) {
1535 _material2min[type][2] = iZ;
1536 }
1537 if (iX > _material2max[type][0]) {
1538 _material2max[type][0] = iX;
1539 }
1540 if (iY > _material2max[type][1]) {
1541 _material2max[type][1] = iY;
1542 }
1543 if (iZ > _material2max[type][2]) {
1544 _material2max[type][2] = iZ;
1545 }
1546 }
1547}
1548
1549// This function compares two discrete normals (discreteNormal, discreteNormal2) in case of a duplicate assignment of boundary types.
1550// The goal of this function is to combine these special boundaryVoxels to an existing one (in this case boundary or externalEdge) according to
1551// the x-, y- and z-values of their discrete normals.
1552// In the following the algorithm is declared only for the x value, but it is also valid for the y and z values.
1553//
1554// for x1 = x2, the new value of x is x1 (1)
1555// for x1*x2 = -1, the new value of x is 0 (2)
1556// for x1*x2 = 0, the new value is 0 (3)
1557//
1558// It may be possible that all three values equal 0. To avoid that the values are tested again (x²+y²+z²==0) and the loosest assumption (3) is
1559// redefined to.
1560//
1561// If x,y and z == 0 --> find those x,y or z which are 0 because of (3) and redefine them to the value !=0
1562//
1563// Additionally the calculated entries are multiplied with (-1) to get the right existing boundary.
1564
1565template<typename T>
1566std::vector<int> BlockGeometryStatistics3D<T>::checkExtraBoundary(
1567 std::vector<int> discreteNormal, std::vector<int> discreteNormal2)
1568{
1569 update();
1570 return const_this->checkExtraBoundary( discreteNormal, discreteNormal2);
1571}
1572
1573template<typename T>
1574std::vector<int> BlockGeometryStatistics3D<T>::checkExtraBoundary(
1575 std::vector<int> discreteNormal, std::vector<int> discreteNormal2) const
1576{
1577 return discreteNormal;
1578}
1579
1580} // namespace olb
1581
1582#endif
blockGeometryStatistics3D.h
Representation of a statistic for a 3D geometry – header file.
blockGeometry.h
Representation of the 2D block geometry view – header file.
olb::BlockGeometryStatistics3D::computeNormal
std::vector< int > computeNormal(int iX, int iY, int iZ)
Returns normal that points into the fluid for paraxial surfaces.
Definition blockGeometryStatistics3D.hh:1316
olb::BlockGeometryStatistics3D::getCenterPhysR
std::vector< T > getCenterPhysR(int material)
Returns the center position.
Definition blockGeometryStatistics3D.hh:260
olb::BlockGeometryStatistics3D::getMaxLatticeR
std::vector< int > getMaxLatticeR(int material)
Returns the max. lattice position in each direction.
Definition blockGeometryStatistics3D.hh:168
olb::BlockGeometryStatistics3D::getNmaterials
int getNmaterials()
Returns the number of different materials.
Definition blockGeometryStatistics3D.hh:89
olb::BlockGeometryStatistics3D::BlockGeometryStatistics3D
BlockGeometryStatistics3D(BlockGeometry< T, 3 > *blockGeometry)
Constructor.
Definition blockGeometryStatistics3D.hh:43
olb::BlockGeometryStatistics3D::getMinPhysR
std::vector< T > getMinPhysR(int material) const
Returns the min. phys position in each direction.
Definition blockGeometryStatistics3D.hh:187
olb::BlockGeometryStatistics3D::getType
std::vector< int > getType(int iX, int iY, int iZ, bool anyNormal=false)
Definition blockGeometryStatistics3D.hh:277
olb::BlockGeometryStatistics3D::getMaxPhysR
std::vector< T > getMaxPhysR(int material) const
Returns the max. phys position in each direction.
Definition blockGeometryStatistics3D.hh:195
olb::BlockGeometryStatistics3D::getNvoxel
std::size_t getNvoxel()
Returns the number of voxels with material!=0.
Definition blockGeometryStatistics3D.hh:133
olb::BlockGeometryStatistics3D::getLatticeExtend
std::vector< T > getLatticeExtend(int material)
Returns the lattice extend as length in each direction.
Definition blockGeometryStatistics3D.hh:203
olb::BlockGeometryStatistics3D::getMaterial2n
std::map< int, std::size_t > getMaterial2n()
Returns the map with the numbers of voxels for each material.
Definition blockGeometryStatistics3D.hh:120
olb::BlockGeometryStatistics3D::getPhysExtend
std::vector< T > getPhysExtend(int material)
Returns the phys extend as length in each direction.
Definition blockGeometryStatistics3D.hh:226
olb::BlockGeometryStatistics3D::getMinLatticeR
std::vector< int > getMinLatticeR(int material)
Returns the min. lattice position in each direction.
Definition blockGeometryStatistics3D.hh:149
olb::BlockGeometryStatistics3D::update
void update(bool verbose=false)
Updates the statistics if it is really needed.
Definition blockGeometryStatistics3D.hh:64
olb::BlockGeometryStatistics3D::computeDiscreteNormal
std::vector< int > computeDiscreteNormal(int material, T maxNorm=1.1)
Returns discrete normal with norm maxNorm that points into the fluid for paraxial surfaces without up...
Definition blockGeometryStatistics3D.hh:1396
olb::BlockGeometryStatistics3D::getStatisticsStatus
bool & getStatisticsStatus()
Read and write access to a flag, which indicates if an uptate is needed (=true)
Definition blockGeometryStatistics3D.hh:51
olb::BlockGeometryStatistics3D::check
bool check(int material, int iX, int iY, int iZ, unsigned offsetX, unsigned offsetY, unsigned offsetZ)
Definition blockGeometryStatistics3D.hh:1429
olb::BlockGeometryStatistics3D::getPhysRadius
std::vector< T > getPhysRadius(int material)
Returns the phys radius as length in each direction.
Definition blockGeometryStatistics3D.hh:243
olb::BlockGeometry
Representation of a block geometry.
Definition blockGeometry.h:56
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
olb::norm
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
Definition scalarVector.h:66
Generated on Tue Mar 5 2024 09:57:50 for OpenLB by doxygen 1.10.0