29#ifndef BLOCK_GEOMETRY_HH
30#define BLOCK_GEOMETRY_HH
43template<
typename T,
unsigned D>
46 _data(this->getNcells()),
47 _communicatable(_data),
51 clout(std::cout, (
"BlockGeometry" + std::to_string(D) +
"D"))
53 _statistics.update(
false);
56template<
typename T,
unsigned D>
62template<
typename T,
unsigned D>
68template<
typename T,
unsigned D>
71 return _cuboid.getOrigin();
74template<
typename T,
unsigned D>
77 return _cuboid.getDeltaR();
80template<
typename T,
unsigned D>
83 return _data[0][this->getCellId(latticeR)];
86template<
typename T,
unsigned D>
89 return _data[0][this->getCellId(latticeR)];
92template<
typename T,
unsigned D>
95 return _data[0][iCell];
98template<
typename T,
unsigned D>
101 if (this->isInside(latticeR)) {
102 return _data[0][this->getCellId(latticeR)];
109template<
typename T,
unsigned D>
113 if (this->isInside(loc)) {
114 return _data[0][this->getCellId(loc)];
121template<
typename T,
unsigned D>
124 set(this->getCellId(latticeR), material);
127template<
typename T,
unsigned D>
130 set(this->getCellId(latticeR), material);
133template<
typename T,
unsigned D>
137 _data[0][iCell] = material;
140template<
typename T,
unsigned D>
144 getPhysR(physR, loc);
148template<
typename T,
unsigned D>
151 _cuboid.getPhysR(physR, loc);
155template<
typename T,
unsigned D>
161template<
typename T,
unsigned D>
164 if constexpr (D == 3) {
165 return Vector<int,3>(this->getNx(), this->getNy(), this->getNz());
169 __builtin_unreachable();
172template<
typename T,
unsigned D>
173template<
typename DESCRIPTOR>
179 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
181 if (get(latticeR) != 0 && (!
util::isContained(bulkMaterials, get(latticeR))) ) {
183 for (
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
184 if (
util::isContained(bulkMaterials, get(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))) ){
195 clout <<
"cleaned "<< counter <<
" outer boundary voxel(s)" << std::endl;
200template<
typename T,
unsigned D>
205 bool toClean =
false;
206 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
209 for(
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
210 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop))) == 0){
221 clout <<
"cleaned "<< counter <<
" outer fluid voxel(s)" << std::endl;
226template<
typename T,
unsigned D>
232 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
233 if (get(latticeR) != 1 && get(latticeR) != 0) {
235 bool var[7] = {
false};
236 for(
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
237 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
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){
252 for(
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
253 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
266 clout <<
"cleaned "<< count2 <<
" inner boundary voxel(s)" << std::endl;
271template<
typename T,
unsigned D>
277 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
278 if (get(latticeR) != 1 && get(latticeR) != 0 && get(latticeR) == fromM) {
280 bool var[7] = {
false};
281 for(
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
282 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
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){
297 for(
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
298 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop)))==1){
311 clout <<
"cleaned "<< count2
312 <<
" inner boundary voxel(s) of Type " << fromM << std::endl;
317template<
typename T,
unsigned D>
320 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
323 getPhysR(physR, latticeR);
324 domain(&output, physR);
331template<
typename T,
unsigned D>
333 std::vector<int> var)
336 for (var[0] = 0; var[0] < this->getNx(); var[0]++) {
337 for (var[1] = 0; var[1] < this->getNy(); var[1]++) {
339 for (var[2] = 0; var[2] < this->getNz(); var[2]++) {
340 found = check(material, var, offset);
343 found = check(material, var, offset);
353template<
typename T,
unsigned D>
355 std::vector<unsigned> offset)
358 for (
int iOffsetX = -offset[0]; iOffsetX <= (int) offset[0]; ++iOffsetX) {
359 for (
int iOffsetY = -offset[1]; iOffsetY <= (int) offset[1]; ++iOffsetY) {
361 for (
int iOffsetZ = -offset[2]; iOffsetZ <= (int) offset[2]; ++iOffsetZ) {
362 if (getMaterial({var[0] + iOffsetX, var[1] + iOffsetY, var[2] + iOffsetZ}) != material) {
367 if (getMaterial({var[0] + iOffsetX, var[1] + iOffsetY}) != material) {
376template<
typename T,
unsigned D>
381 bool errorFound =
false;
382 this->forSpatialLocations([&](
LatticeR<D> latticeR)
384 if (get(latticeR) == 0) {
386 for(
int iPop = 1; iPop < DESCRIPTOR::q; iPop++){
387 if(getMaterial(latticeR + (descriptors::c<DESCRIPTOR>(iPop))) == 1){
399 clout <<
"error!" << std::endl;
402 clout <<
"the model is correct!" << std::endl;
408template<
typename T,
unsigned D>
411 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
412 if (get(latticeR) == fromM) {
418template<
typename T,
unsigned D>
422 this->forSpatialLocations([&](
LatticeR<D> latticeR) {
423 if (get(latticeR) == fromM) {
424 getPhysR(physR, latticeR);
426 condition(inside, physR);
434template<
typename T,
unsigned D>
437 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
438 if (get(latticeR) == fromM) {
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) {
451 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY}) != fromM) {
452 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY}) != 1245) {
467template<
typename T,
unsigned D>
469 std::vector<int> testDirection)
471 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
472 if (get(latticeR) == fromM) {
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) {
486 if (iOffsetX!=0 || iOffsetY!=0) {
487 if (getMaterial({latticeR[0] + iOffsetX, latticeR[1] + iOffsetY}) != testM) {
502template<
typename T,
unsigned D>
506 rename(fromM, toM, condition);
509 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
510 if (get(latticeR) == toM) {
511 getPhysR(physR, latticeR);
513 condition(inside, physR);
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);
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);
533template<
typename T,
unsigned D>
537 rename(fromM, toM, condition);
538 std::vector<int> testDirection = getStatistics().computeDiscreteNormal(toM);
541 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
542 if (get(latticeR) == toM) {
543 getPhysR(physR, latticeR);
545 condition(inside, physR);
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);
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);
565template<
typename T,
unsigned D>
569 this->forCoreSpatialLocations([&](
LatticeR<D> latticeR) {
570 getPhysR(physR, latticeR);
572 condition(inside, physR);
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()){
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));
582 set({latticeR[0] + i * discreteNormal[0], latticeR[1] + i * discreteNormal[1]}, get(latticeR));
590template<
typename T,
unsigned D>
592 std::vector<int> offset, std::map<std::vector<int>,
int>* tmp)
594 std::map<std::vector<int>,
int> tmp2;
595 bool firstCall =
false;
596 if (tmp ==
nullptr) {
601 if (getMaterial(seed) == fromM) {
602 std::vector<int> found;
603 found.push_back(seed[0]);
604 found.push_back(seed[1]);
606 found.push_back(seed[2]);
607 if (tmp->count(found) == 0) {
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);
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);
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);
623 if (tmp->count(found) == 0) {
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);
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);
637 std::map<std::vector<int>,
int>::iterator iter;
638 for (iter = tmp->begin(); iter != tmp->end(); iter++) {
640 set((iter->first)[0],(iter->first)[1],(iter->first)[2], toM);
642 set((iter->first)[0],(iter->first)[1], toM);
649template<
typename T,
unsigned D>
652 for (
int x = min[0]; x <= max[0]; x++) {
656 for (
int y = min[1]; y <= max[1]; y++) {
658 for (
int z = min[2]; z <= max[2]; z++) {
659 clout << getMaterial({x, y, z}) <<
" ";
661 if (max[1] - min[1] != 0 && max[2] - min[2] != 0) {
665 clout << getMaterial({x, y}) <<
" ";
668 if (max[0] - min[0] != 0) {
675template<
typename T,
unsigned D>
678 assert(direction >= 0 && direction <= 2);
682 printLayer({layer, 0, 0},{layer, this->getNy() - 1, this->getNz() - 1}, linenumber);
685 printLayer({0, layer, 0}, {this->getNx() - 1, layer, this->getNz() - 1}, linenumber);
688 printLayer({0, 0, layer}, {this->getNx() - 1, this->getNy() - 1, layer}, linenumber);
694 printLayer({layer, 0}, {layer, this->getNy() - 1}, linenumber);
697 printLayer({0, layer}, {this->getNx() - 1, layer}, linenumber);
703template<
typename T,
unsigned D>
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++) {
710 for (
int z = loc[2] - 1; z <= loc[2] + 1; z++) {
711 clout << getMaterial({x, y, z}) <<
" ";
715 clout << getMaterial({x, y}) <<
" ";
723template<
typename T,
unsigned D>
726 _statistics.getStatisticsStatus() =
true;
729template<
typename T,
unsigned D>
735template<
typename T,
unsigned D>
738 return _data.getSerializableSize();
741template<
typename T,
unsigned D>
744 std::size_t currentBlock = 0;
745 bool* dataPtr =
nullptr;
747 this->registerSerializableOfConstSize(iBlock, sizeBlock, currentBlock, dataPtr, _data, loadingMode);
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.
IndicatorF3D is an application from .
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
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)
Top level namespace for all of OpenLB.
std::conditional_t< D==2, BlockGeometryStatistics2D< T >, BlockGeometryStatistics3D< T > > BlockGeometryStatistics
std::conditional_t< D==2, IndicatorF2D< T >, IndicatorF3D< T > > IndicatorF
std::conditional_t< D==2, Cuboid2D< T >, Cuboid3D< T > > Cuboid