diff --git a/apps/lbmMultiRes/flowOverShape.h b/apps/lbmMultiRes/flowOverShape.h index aecd70e6..ec843d07 100644 --- a/apps/lbmMultiRes/flowOverShape.h +++ b/apps/lbmMultiRes/flowOverShape.h @@ -5,14 +5,15 @@ #include "init.h" -#include "igl/AABB.h" -#include "igl/WindingNumberAABB.h" + #include "igl/max.h" #include "igl/min.h" #include "igl/read_triangle_mesh.h" #include "igl/remove_unreferenced.h" #include "igl/writeOBJ.h" +#include "igl/AABB.h" + template void initFlowOverShape(Neon::domain::mGrid& grid, Neon::domain::mGrid::Field& sumStore, @@ -278,15 +279,6 @@ void flowOverSphere(const Neon::Backend backend, } -template -inline auto winding_number_aabbc(const Eigen::MatrixBase& V, - const Eigen::MatrixBase& F) -{ - return igl::WindingNumberAABB< - Eigen::Matrix, - DerivedV, - DerivedF>(V, F); -} template void flowOverMesh(const Neon::Backend backend, const Params& params) @@ -294,11 +286,14 @@ void flowOverMesh(const Neon::Backend backend, static_assert(std::is_same_v || std::is_same_v); //define the gird and the box that will encompass the mesh - Neon::index_3d gridDim(19 * params.scale, 10 * params.scale, 10 * params.scale); - - Eigen::RowVector3d meshBoxDim(4 * params.scale, 4 * params.scale, 4 * params.scale); + Neon::index_3d gridDim(19 * params.scale, 10 * params.scale, 10 * params.scale); + Eigen::RowVector3d meshBoxDim(3.5 * params.scale, 3.5 * params.scale, 3.5 * params.scale); Eigen::RowVector3d meshBoxCenter(5 * params.scale, 5 * params.scale, 5 * params.scale); + //Neon::index_3d gridDim(1.5 * 29 * params.scale, 10 * params.scale, 10 * params.scale); + //Eigen::RowVector3d meshBoxDim(2 * params.scale, 2 * params.scale, 2 * params.scale); + //Eigen::RowVector3d meshBoxCenter(20 * params.scale, 10, 5 * params.scale); + //read the mesh and scale it such that it fits inside meshBox Eigen::MatrixXi faces; @@ -309,18 +304,50 @@ void flowOverMesh(const Neon::Backend backend, Eigen::VectorXi _1, _2; igl::remove_unreferenced(Eigen::MatrixXd(vertices), Eigen::MatrixXi(faces), vertices, faces, _1, _2); + //mesh bounding box using the mesh coordinates Eigen::RowVectorXd bbMax, bbMin; Eigen::RowVectorXi bbMaxI, bbMinI; igl::max(vertices, 1, bbMax, bbMaxI); igl::min(vertices, 1, bbMin, bbMinI); - //translate and scale the mesh + //center the mesh vertices.rowwise() -= ((bbMin + bbMax) / 2.0); + + //scale double scaling_factor = meshBoxDim.maxCoeff() / (bbMax - bbMin).maxCoeff(); vertices *= scaling_factor; + + auto toRad = [](double t) { return t * (3.14159265358979311600 / 180); }; + + Eigen::Matrix3d rotationX; + //rotate about x-axis by thetaX degrees + rotationX << 1, 0, 0, + 0, std::cos(toRad(params.thetaX)), -std::sin(toRad(params.thetaX)), + 0, std::sin(toRad(params.thetaX)), std::cos(toRad(params.thetaX)); + + //rotate about y-axis by thetaY degrees + Eigen::Matrix3d rotationY; + rotationY << std::cos(toRad(params.thetaY)), 0, std::sin(toRad(params.thetaY)), + 0, 1, 0, + -std::sin(toRad(params.thetaY)), 0, std::cos(toRad(params.thetaY)); + + //rotate about z-axis by thetaZ degrees + Eigen::Matrix3d rotationZ; + rotationZ << std::cos(toRad(params.thetaZ)), -std::sin(toRad(params.thetaZ)), 0, + std::sin(toRad(params.thetaZ)), std::cos(toRad(params.thetaZ)), 0, + 0, 0, 1; + + Eigen::Matrix3d rotation = rotationX * rotationY * rotationZ; + vertices = vertices * rotation; + + //translate vertices.rowwise() += meshBoxCenter; + //calc the bounding box again + igl::max(vertices, 1, bbMax, bbMaxI); + igl::min(vertices, 1, bbMin, bbMinI); + igl::writeOBJ("scaled.obj", vertices, faces); @@ -328,41 +355,88 @@ void flowOverMesh(const Neon::Backend backend, polyscopeAddMesh(params.meshFile, faces, vertices); #endif - NEON_INFO("Winding Number init"); - //get the mesh ready for query (inside/outside) using libigl winding number - auto wnAABB = winding_number_aabbc(vertices, faces); - wnAABB.set_method(igl::WindingNumberMethod::APPROX_SIMPLE_WINDING_NUMBER_METHOD); - wnAABB.grow(); + NEON_INFO("AABB init"); + + igl::AABB aabb; + aabb.init(vertices, faces); + //define the grid - int depth = 3; + int depth = 3; + + auto distToMesh = [&](const Neon::index_3d idx) { + Eigen::Matrix p(idx.x, idx.y, idx.z), c; + int face_index(0); + aabb.squared_distance(vertices, faces, p, face_index, c); + return (p - c).norm(); + }; + + std::vector> + activeCellLambda = + {[&](const Neon::index_3d idx) -> bool { + return idx.x >= 2 * params.scale && idx.x < 8 * params.scale && + idx.y >= 3 * params.scale && idx.y < 7 * params.scale && + idx.z >= 3 * params.scale && idx.z < 7 * params.scale; + //return distToMesh(idx) < params.scale; + }, + [&](const Neon::index_3d idx) -> bool { + return idx.x >= params.scale && idx.x < 13 * params.scale && + idx.y >= 2 * params.scale && idx.y < 8 * params.scale && + idx.z >= 2 * params.scale && idx.z < 8 * params.scale; + //return distToMesh(idx) < 2 * params.scale; + }, + [&](const Neon::index_3d idx) -> bool { + return true; + }}; + + //int depth = 5; + // + //std::vector> + // activeCellLambda = + // {[&](const Neon::index_3d idx) -> bool { + // return idx.x >= (0.1 * 29 * params.scale) + 4.3 * params.scale && idx.x < (0.1 * 29 * params.scale) + 10 * params.scale && + // idx.y < 1.2 * params.scale && + // idx.z >= 3 * params.scale && idx.z < 7 * params.scale; + // }, + // [&](const Neon::index_3d idx) -> bool { + // return idx.x >= (0.1 * 29 * params.scale) + 3.0 * params.scale && idx.x < (0.1 * 29 * params.scale) + 14 * params.scale && + // idx.y < 3 * params.scale && + // idx.z >= 3 * params.scale && idx.z < 7 * params.scale; + // }, + // [&](const Neon::index_3d idx) -> bool { + // return idx.x >= (0.1 * 29 * params.scale) + 2.2 * params.scale && idx.x < (0.1 * 29 * params.scale) + 19.5 * params.scale && + // idx.y < 4.3 * params.scale && + // idx.z >= 3 * params.scale && idx.z < 7 * params.scale; + // }, + // [&](const Neon::index_3d idx) -> bool { + // return idx.x >= (0.1 * 29 * params.scale) + 0.7 * params.scale && idx.x < (0.1 * 29 * params.scale) + 25 * params.scale && + // idx.y < 6.2 * params.scale && + // idx.z >= 2 * params.scale && idx.z < 8 * params.scale; + // }, + // [&](const Neon::index_3d idx) -> bool { + // return true; + // }}; + const Neon::mGridDescriptor<1> descriptor(depth); + NEON_INFO("Create mGrid"); + Neon::domain::mGrid grid( - backend, gridDim, - {[&](const Neon::index_3d idx) -> bool { - return idx.x >= 2 * params.scale && idx.x < 8 * params.scale && - idx.y >= 3 * params.scale && idx.y < 7 * params.scale && - idx.z >= 3 * params.scale && idx.z < 7 * params.scale; - }, - [&](const Neon::index_3d idx) -> bool { - return idx.x >= params.scale && idx.x < 13 * params.scale && - idx.y >= 2 * params.scale && idx.y < 8 * params.scale && - idx.z >= 2 * params.scale && idx.z < 8 * params.scale; - }, - [&](const Neon::index_3d idx) -> bool { - return true; - }}, - Neon::domain::Stencil::s19_t(false), descriptor); + backend, gridDim, activeCellLambda, Neon::domain::Stencil::s19_t(false), descriptor); //LBM problem - const T uin = 0.04; - const T clength = T((meshBoxDim.minCoeff() / 2) / (1 << (depth - 1))); + const T uin = 0.04; + //const T clength = T((meshBoxDim.minCoeff() / 2) / (1 << (depth - 1))); + + //the plane wing extends along the z direction + const T span = (bbMax - bbMin).z(); + const T clength = (span / 2.f) / (1 << (depth - 1)); const T visclb = uin * clength / static_cast(params.Re); const T omega = 1.0 / (3. * visclb + 0.5); const Neon::double_3d inletVelocity(uin, 0., 0.); + //NEON_INFO("Writing Test file"); //auto test = grid.newField("test", 1, 0); //test.ioToVtk("Test", true, true, true, true, {-1, -1, 1}); //exit(0); @@ -383,9 +457,19 @@ void flowOverMesh(const Neon::Backend backend, voxelGlobalLocation.z < meshBoxCenter.z() - meshBoxDim.z() / 2 || voxelGlobalLocation.z >= meshBoxCenter.z() + meshBoxDim.z() / 2) { in(cell, 0) = 0; } else { - const double voxelSpacing = 0.5 * double(grid.getDescriptor().getSpacing(l - 1)); - Eigen::RowVector3d point(voxelGlobalLocation.x + voxelSpacing, voxelGlobalLocation.y + voxelSpacing, voxelGlobalLocation.z + voxelSpacing); - in(cell, 0) = int8_t(wnAABB.winding_number(point) + std::numeric_limits::epsilon() > 1); + const double voxelSpacing = 0.5 * double(grid.getDescriptor().getSpacing(l - 1)); + + const double centerToCornerDistSqr = (3.0 / 4.0) * (2.0 * voxelSpacing); + + Eigen::Matrix p(voxelGlobalLocation.x + voxelSpacing, voxelGlobalLocation.y + voxelSpacing, voxelGlobalLocation.z + voxelSpacing), c; + + int face_index(0); + + aabb.squared_distance(vertices, faces, p, face_index, c); + + double sqr_dist = (p - c).norm(); + + in(cell, 0) = int8_t(sqr_dist + std::numeric_limits::epsilon() < centerToCornerDistSqr); } } else { in(cell, 0) = 0; diff --git a/apps/lbmMultiRes/lbmMultiRes.cu b/apps/lbmMultiRes/lbmMultiRes.cu index 66bd37f7..af1953d8 100644 --- a/apps/lbmMultiRes/lbmMultiRes.cu +++ b/apps/lbmMultiRes/lbmMultiRes.cu @@ -1,7 +1,7 @@ #include "Neon/core/tools/clipp.h" -#define BGK -//#define KBC +//#define BGK +#define KBC #include "Neon/Neon.h" #include "Neon/Report.h" @@ -30,7 +30,11 @@ struct Params int sliceX = -1; int sliceY = -1; int sliceZ = 1; + double thetaX = 0; + double thetaY = 0; + double thetaZ = 0; bool vtk = false; + bool binary = false; bool gui = false; int scale = 2; std::string dataType = "float"; @@ -61,10 +65,15 @@ int main(int argc, char** argv) clipp::option("--sliceY") & clipp::integer("sliceY", params.sliceY) % "Slice along Y for output images/VTK", clipp::option("--sliceZ") & clipp::integer("sliceZ", params.sliceZ) % "Slice along Z for output images/VTK", + clipp::option("--thetaX") & clipp::value("thetaX", params.thetaX) % "Angle (degree) of rotation of the input model along X axis", + clipp::option("--thetaY") & clipp::value("thetaY", params.thetaY) % "Angle (degree) of rotation of the input model along Y axis", + clipp::option("--thetaZ") & clipp::value("thetaZ", params.thetaZ) % "Angle (degree) of rotation of the input model along Z axis", + ((clipp::option("--benchmark").set(params.benchmark, true) % "Run benchmark mode") | (clipp::option("--visual").set(params.benchmark, false) % "Run export partial data")), clipp::option("--vtk").set(params.vtk, true) % "Output VTK files. Active only with if 'visual' is true", + clipp::option("--binary").set(params.binary, true) % "Output binary (down-sampled) files. Active only with if 'visual' is true", clipp::option("--gui").set(params.gui, true) % "Show Polyscope gui. Active only with if 'visual' is true", clipp::option("--freq") & clipp::integers("freq", params.freq) % "Output frequency (only works with visual mode)", diff --git a/apps/lbmMultiRes/lbmMultiRes.h b/apps/lbmMultiRes/lbmMultiRes.h index 2b013d95..1143f422 100644 --- a/apps/lbmMultiRes/lbmMultiRes.h +++ b/apps/lbmMultiRes/lbmMultiRes.h @@ -369,7 +369,7 @@ void runNonUniformLBM(Neon::domain::mGrid& grid, suffix << std::setw(precision) << std::setfill('0') << t; std::string fileName = "Velocity_" + suffix.str(); - postProcess(grid, depth, fout, cellType, vel, rho, slice, fileName, params.vtk && t != 0, psDrawable, psHex, psHexVert); + postProcess(grid, depth, fout, cellType, vel, rho, slice, fileName, params.vtk && t != 0, params.binary && t != 0, psDrawable, psHex, psHexVert); #ifdef NEON_USE_POLYSCOPE postProcessPolyscope(psDrawable, vel, psColor, fileName, params.gui, t == 0); #endif @@ -486,7 +486,7 @@ void runNonUniformLBM(Neon::domain::mGrid& grid, report.addMember("ENumVoxels", gridDim.rMul()); //output - report.write("MultiResLBM_" + reportSuffix(), true); + //report.write("MultiResLBM_" + reportSuffix(), true); //post process if (!params.benchmark) { @@ -494,12 +494,12 @@ void runNonUniformLBM(Neon::domain::mGrid& grid, std::ostringstream suffix; suffix << std::setw(precision) << std::setfill('0') << params.numIter; std::string fileName = "Velocity_" + suffix.str(); - postProcess(grid, depth, fout, cellType, vel, rho, slice, fileName, params.vtk, psDrawable, psHex, psHexVert); + postProcess(grid, depth, fout, cellType, vel, rho, slice, fileName, params.vtk, params.binary, psDrawable, psHex, psHexVert); #ifdef NEON_USE_POLYSCOPE postProcessPolyscope(psDrawable, vel, psColor, fileName, params.gui, false); #endif } else if (verify) { - postProcess(grid, depth, fout, cellType, vel, rho, slice, "", false, psDrawable, psHex, psHexVert); + postProcess(grid, depth, fout, cellType, vel, rho, slice, "", false, false, psDrawable, psHex, psHexVert); } if (verify) { @@ -510,60 +510,4 @@ void runNonUniformLBM(Neon::domain::mGrid& grid, params.numIter, velocity.x); } - - if (!params.benchmark) { - NEON_INFO("Started Binary VTK"); - - //the level at which we will do the sampling - const int theLevel = depth - 1; - - const Neon::index_4d grid4D(gridDim.x / (1 << theLevel), gridDim.y / (1 << theLevel), gridDim.z / (1 << theLevel), 3); - std::vector ioBuffer(grid4D.x * grid4D.y * grid4D.z * grid4D.w); - float* ioBufferPtr = ioBuffer.data(); - - for (int l = 0; l < grid.getDescriptor().getDepth(); ++l) { - - grid.newContainer("Viz", l, [=](Neon::set::Loader& loader) { - auto& v = vel.load(loader, l, Neon::MultiResCompute::MAP); - - return [=](const typename Neon::domain::mGrid::Idx& cell) mutable { - if (!v.hasChildren(cell)) { - Neon::index_3d loc = v.getGlobalIndex(cell); - - if (loc.x % (1 << theLevel) != 0 || loc.y % (1 << theLevel) != 0 || loc.z % (1 << theLevel) != 0) { - return; - } - - loc.x /= (1 << theLevel); - loc.y /= (1 << theLevel); - loc.z /= (1 << theLevel); - - - for (int c = 0; c < 3; ++c) { - - const float val = v(cell, c); - - - Neon::index_4d loc4D(loc.x, loc.y, loc.z, c); - ioBufferPtr[loc4D.mPitch(grid4D)] = val; - } - } - }; - }) - .run(0); - } - - { - Neon::index_3d nodeDim(grid4D.x + 1, grid4D.y + 1, grid4D.z + 1); - Neon::IoToVTK io("BinVelocity", nodeDim, {1, 1, 1}, {0, 0, 0}, Neon::IoFileType::BINARY); - - io.addField([=](Neon::index_3d idx, int card) { - Neon::index_4d loc4D(idx.x, idx.y, idx.z, card); - return ioBufferPtr[loc4D.mPitch(grid4D)]; - }, - 3, "V", Neon::ioToVTKns::VtiDataType_e::voxel); - } - - NEON_INFO("Finished Binary VTK"); - } } diff --git a/apps/lbmMultiRes/postProcess.h b/apps/lbmMultiRes/postProcess.h index 2526141e..de7dac06 100644 --- a/apps/lbmMultiRes/postProcess.h +++ b/apps/lbmMultiRes/postProcess.h @@ -22,6 +22,7 @@ void postProcess(Neon::domain::mGrid& const Neon::int8_3d slice, std::string fileName, bool outputFile, + bool outoutFileBinary, const std::vector>& psDrawable, const std::vector>& psHex, const std::vector& psHexVert) @@ -41,7 +42,7 @@ void postProcess(Neon::domain::mGrid& return [=] NEON_CUDA_HOST_DEVICE(const typename Neon::domain::mGrid::Idx& cell) mutable { if (!pop.hasChildren(cell)) { - if (type(cell, 0) == CellType::bulk ) { + if (type(cell, 0) == CellType::bulk) { //fin T ins[Q]; @@ -144,10 +145,10 @@ void postProcess(Neon::domain::mGrid& for (size_t t = 0; t < psDrawable.size(); ++t) { const auto id = psDrawable[t].first; int level = psDrawable[t].second; - + for (int d = 0; d < 3; ++d) { T v = vel(id, d, level); - file << v << " "; + file << v << " "; } file << "\n"; @@ -159,8 +160,71 @@ void postProcess(Neon::domain::mGrid& //file << c << "\n"; } + + file << "SCALARS Level float 1 \n"; + file << "LOOKUP_TABLE default \n"; + for (size_t t = 0; t < psDrawable.size(); ++t) { + int level = psDrawable[t].second; + file << level << "\n"; + } + file.close(); } + + if (outoutFileBinary) { + //the level at which we will do the sampling + const int depth = grid.getDescriptor().getDepth(); + const auto gridDim = grid.getDimension(); + const int theLevel = depth - 1; + + const Neon::index_4d grid4D(gridDim.x / (1 << theLevel), gridDim.y / (1 << theLevel), gridDim.z / (1 << theLevel), 3); + std::vector ioBuffer(grid4D.x * grid4D.y * grid4D.z * grid4D.w); + float* ioBufferPtr = ioBuffer.data(); + + for (int l = 0; l < grid.getDescriptor().getDepth(); ++l) { + + grid.newContainer("Viz", l, [=](Neon::set::Loader& loader) { + auto& v = vel.load(loader, l, Neon::MultiResCompute::MAP); + + return [=](const typename Neon::domain::mGrid::Idx& cell) mutable { + if (!v.hasChildren(cell)) { + Neon::index_3d loc = v.getGlobalIndex(cell); + + if (loc.x % (1 << theLevel) != 0 || loc.y % (1 << theLevel) != 0 || loc.z % (1 << theLevel) != 0) { + return; + } + + loc.x /= (1 << theLevel); + loc.y /= (1 << theLevel); + loc.z /= (1 << theLevel); + + + for (int c = 0; c < 3; ++c) { + + const float val = v(cell, c); + + + Neon::index_4d loc4D(loc.x, loc.y, loc.z, c); + ioBufferPtr[loc4D.mPitch(grid4D)] = val; + } + } + }; + }) + .run(0); + } + + { + Neon::index_3d nodeDim(grid4D.x + 1, grid4D.y + 1, grid4D.z + 1); + + Neon::IoToVTK io("Bin" + fileName, nodeDim, {1, 1, 1}, {0, 0, 0}, Neon::IoFileType::BINARY); + + io.addField([=](Neon::index_3d idx, int card) { + Neon::index_4d loc4D(idx.x, idx.y, idx.z, card); + return ioBufferPtr[loc4D.mPitch(grid4D)]; + }, + 3, "V", Neon::ioToVTKns::VtiDataType_e::voxel); + } + } } diff --git a/apps/lbmMultiRes/scripts/MultiResNeon_vs_ghia1982.png b/apps/lbmMultiRes/scripts/MultiResNeon_vs_ghia1982.png index 431ceb45..9c04c98b 100644 Binary files a/apps/lbmMultiRes/scripts/MultiResNeon_vs_ghia1982.png and b/apps/lbmMultiRes/scripts/MultiResNeon_vs_ghia1982.png differ diff --git a/apps/lbmMultiRes/scripts/UniformNeon_vs_ghia1982.png b/apps/lbmMultiRes/scripts/UniformNeon_vs_ghia1982.png index af27920a..a8e02af4 100644 Binary files a/apps/lbmMultiRes/scripts/UniformNeon_vs_ghia1982.png and b/apps/lbmMultiRes/scripts/UniformNeon_vs_ghia1982.png differ diff --git a/apps/lbmMultiRes/scripts/plot.py b/apps/lbmMultiRes/scripts/plot.py index ad068db5..64738507 100644 --- a/apps/lbmMultiRes/scripts/plot.py +++ b/apps/lbmMultiRes/scripts/plot.py @@ -1,45 +1,65 @@ +#https://help.altair.com/hwcfdsolvers/nfx/topics/nanofluidx/lid_driven_cavity_2d_r.htm import os import numpy as np import matplotlib.pyplot as plt +from matplotlib import rcParams +rcParams['font.family'] = ['Palatino Linotype', 'serif'] +rcParams['font.size'] = 10 +rcParams["mathtext.default"] = 'regular' def plotme(dat_filename_X="", dat_filename_Y="", label="", png_filename="", scale=1): - fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(4, 3), dpi=200) - if (dat_filename_Y!=""): - - - y_neon, u_neon = np.loadtxt(dat_filename_Y, unpack=True, usecols=(0, 1)) - y_neon -=0.5 - axes.plot(y_neon, scale*u_neon, 'y.-', label=label) - - - y_ref, u_ref = np.loadtxt('ghia1982.dat', unpack=True, skiprows=2, usecols=(0, 1)) - y_ref -=0.5 - axes.plot(y_ref, u_ref, 'bs', markerfacecolor='none', label='Ghia et al. 1982') - - plt.xticks([-0.5, -0.25, 0.0, 0.25, 0.5]) + fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(4, 3), dpi=1000) if (dat_filename_X!=""): x_neon, v_neon = np.loadtxt(dat_filename_X, unpack=True, usecols=(0, 1)) x_neon -= 0.5 - axes.plot(x_neon, scale*v_neon, 'g.-', label=label) + x_neon *= 2.0 + ax1.plot(x_neon, scale*v_neon, 'k-', label=label) x_ref, v_ref = np.loadtxt('ghia1982.dat', unpack=True, skiprows=2, usecols=(6, 7)) x_ref -= 0.5 - axes.plot(x_ref, v_ref, 'ks', markerfacecolor='none', - label='Ghia et al. 1982') + x_ref *= 2.0 + ax1.plot(x_ref, v_ref, 'o', markerfacecolor='none', + label='Ghia 1982') + + ax1.set_xlabel(r'x') + ax1.set_ylabel('v/u$_{lid}$') + ax1.legend(loc = "upper left") + + if (dat_filename_Y!=""): + ax2 = ax1 .twinx() + ax3 = ax1 .twiny() + + y_neon, u_neon = np.loadtxt(dat_filename_Y, unpack=True, usecols=(0, 1)) + y_neon -=0.5 + y_neon *= 2.0 + ax2.plot(scale*u_neon, y_neon, 'k-', label=label) + + ax3.set_xlim(ax2.get_ylim()) + + y_ref, u_ref = np.loadtxt('ghia1982.dat', unpack=True, skiprows=2, usecols=(0, 1)) + y_ref -=0.5 + y_ref *= 2.0 + ax3.plot(u_ref, y_ref, 'o', markerfacecolor='none', label='Ghia 1982') + + ax2.set_ylabel('u/u$_{lid}$') + ax3.set_xlabel(r'y') + + #h1, l1 = ax1.get_legend_handles_labels() + #h2, l2 = ax2.get_legend_handles_labels() + #h3, l3 = ax3.get_legend_handles_labels() + #ax1.legend(h1+h2+h3, l1+l2+l3, loc=2) - axes.legend() - axes.set_xlabel(r'Distance') - axes.set_ylabel(r'Velocity') + plt.tight_layout() plt.savefig(png_filename) plotme(dat_filename_X='NeonMultiResLBM_5000_X.dat', dat_filename_Y='NeonMultiResLBM_5000_Y.dat', - label='Grid Refinement LBM', + label='Ours', png_filename='MultiResNeon_vs_ghia1982.png', scale=1) diff --git a/libNeonDomain/src/domain/details/mGrid/mGrid.cpp b/libNeonDomain/src/domain/details/mGrid/mGrid.cpp index ae200ae5..aaf57ab8 100644 --- a/libNeonDomain/src/domain/details/mGrid/mGrid.cpp +++ b/libNeonDomain/src/domain/details/mGrid/mGrid.cpp @@ -122,6 +122,21 @@ mGrid::mGrid( } + if (containVoxels) { + for (int z = 0; z < refFactor; z++) { + for (int y = 0; y < refFactor; y++) { + for (int x = 0; x < refFactor; x++) { + + const Neon::int32_3d voxel = mData->mDescriptor.parentToChild(blockOrigin, l, {x, y, z}); + + if (voxel < domainSize) { + setLevelBitMask(l, {bx, by, bz}, {x, y, z}); + } + } + } + } + } + if (containVoxels) { //if the block contains voxels, it should activate itself //find its corresponding index within the next level