diff --git a/Source/Constants.H b/Source/Constants.H index e0a4aaae..8ce0a5ed 100644 --- a/Source/Constants.H +++ b/Source/Constants.H @@ -46,7 +46,7 @@ static constexpr int D_Q_CORR_X_IDX = 4; static constexpr int D_Q_CORR_Y_IDX = 5; static constexpr int D_Q_CORR_Z_IDX = 6; -static constexpr int N_IS_FLUID = 2; +static constexpr int N_IS_FLUID = 3; static constexpr amrex::Real ROOT3 = static_cast(1.732050807568877293527446341); diff --git a/Source/EB.H b/Source/EB.H index c34939a9..636599ec 100644 --- a/Source/EB.H +++ b/Source/EB.H @@ -29,6 +29,11 @@ namespace lbm { void initialize_eb(const amrex::Geometry& geom, const int max_level); void initialize_from_stl( const amrex::Geometry& geom, amrex::iMultiFab& is_fluid); +std::vector +read_crack_file(const std::string& filename, int nx, int ny, int nz); + +void generate_voxel_cracks( + const amrex::Geometry& geom, amrex::iMultiFab& is_fluid); } // namespace lbm #endif diff --git a/Source/EB.cpp b/Source/EB.cpp index e951f018..e888a25b 100644 --- a/Source/EB.cpp +++ b/Source/EB.cpp @@ -1,5 +1,7 @@ #include "EB.H" #include "Geometry.H" +#include // for std::memcpy +#include // for CSV parsing namespace lbm { void initialize_eb(const amrex::Geometry& geom, const int max_level) @@ -27,12 +29,15 @@ void initialize_eb(const amrex::Geometry& geom, const int max_level) amrex::Vector amrex_defaults( {"all_regular", "box", "cylinder", "plane", "sphere", "torus", "parser", "stl"}); - if (!(std::find(amrex_defaults.begin(), amrex_defaults.end(), geom_type) != - amrex_defaults.end())) { + if (std::find(amrex_defaults.begin(), amrex_defaults.end(), geom_type) == + amrex_defaults.end()) { std::unique_ptr geometry( lbm::Geometry::create(geom_type)); geometry->build(geom, max_coarsening_level); } else { + // For all AMReX default types (including voxel_cracks), use standard + // build voxel_cracks will override the m_is_fluid in + // initialize_from_stl amrex::EB2::Build(geom, max_level, max_level); } } @@ -84,11 +89,210 @@ void initialize_from_stl( static_cast(marker_arrs[nbx](i, j, k, 0)); }); amrex::Gpu::synchronize(); - } else if ((!name.empty()) && (geom_type != "all_regular")) { + } + + // Check for voxel crack generation flag + int use_voxel_cracks = 0; + pp.query("use_voxel_cracks", use_voxel_cracks); + if (use_voxel_cracks != 0) { + amrex::Print() << "Using voxel crack generation" << std::endl; + generate_voxel_cracks(geom, is_fluid); + return; + } + + if ((!name.empty()) && (geom_type != "all_regular")) { amrex::Abort( "LBM::initialize_from_stl() geom_type should be all_regular to " "avoid issues"); } } +std::vector +read_crack_file(const std::string& filename, int nx, int ny, int nz) +{ + BL_PROFILE("LBM::read_crack_file()"); + + // Auto-detect file format + bool is_csv = + (filename.size() >= 4 && + filename.substr(filename.size() - 4) == ".csv"); + + // Use AMReX's cross-platform file reading utilities + amrex::Vector file_char_ptr; + + // Use AMReX's parallel-safe file reading + try { + amrex::ParallelDescriptor::ReadAndBcastFile(filename, file_char_ptr); + } catch (const std::exception& e) { + amrex::Abort( + "Error reading crack file: " + filename + " - " + + std::string(e.what())); + } + + std::vector crack_data( + static_cast(nx) * static_cast(ny) * + static_cast(nz)); + + if (is_csv) { + // Parse CSV format + if (amrex::ParallelDescriptor::IOProcessor()) { + std::fill( + crack_data.begin(), crack_data.end(), 1); // Initialize as solid + + std::string file_content( + file_char_ptr.data(), file_char_ptr.size()); + std::istringstream iss(file_content); + std::string line; + + // Skip header line + if (!std::getline(iss, line)) { + amrex::Abort("CSV file is empty or corrupted: " + filename); + } + + size_t line_count = 0; + while (std::getline(iss, line)) { + // Skip empty lines + if (line.empty()) { + continue; + } + + std::istringstream line_stream(line); + std::string token; + + // Parse tokens: X, Y, Z (discard), then tag + int tag; + try { + // discard X + if (!std::getline(line_stream, token, ',')) { + continue; + } + // discard Y + if (!std::getline(line_stream, token, ',')) { + continue; + } + // discard Z + if (!std::getline(line_stream, token, ',')) { + continue; + } + // read tag + if (std::getline(line_stream, token, ',')) { + tag = std::stoi(token); + } else { + continue; + } + } catch (const std::exception& e) { + // Skip malformed lines + continue; + } + + // CSV and binary are written in identical sequence by + // mainCrackGenerator.C Both loop: k(Z) -> j(Y) -> i(X), so just + // read sequentially The X,Y,Z coordinates are metadata - what + // matters is the order + + if (line_count < static_cast(nx) * + static_cast(ny) * + static_cast(nz)) { + crack_data[line_count] = static_cast(tag); + } + line_count++; + } + + amrex::Print() << "Successfully read CSV crack file: " << filename + << " (" << line_count << " data points)" + << std::endl; + } + // Broadcast the parsed data to all processors + amrex::ParallelDescriptor::Bcast( + crack_data.data(), crack_data.size() * sizeof(uint16_t), 0); + + } else { + // Parse binary format + size_t expected_size = static_cast(nx) * + static_cast(ny) * + static_cast(nz) * sizeof(uint16_t); + + // ReadAndBcastFile may add extra bytes, so ensure we only use what we + // need + if (static_cast(file_char_ptr.size()) < expected_size) { + amrex::Abort( + "Binary file too small after reading. Expected: " + + std::to_string(expected_size) + " bytes, got: " + + std::to_string(file_char_ptr.size()) + " bytes"); + } + + // Convert char data to uint16_t array + std::memcpy(crack_data.data(), file_char_ptr.data(), expected_size); + + if (amrex::ParallelDescriptor::IOProcessor()) { + amrex::Print() << "Successfully read binary crack file: " + << filename << " (" << expected_size << " bytes)" + << std::endl; + } + } + + return crack_data; +} + +void generate_voxel_cracks( + const amrex::Geometry& geom, amrex::iMultiFab& is_fluid) +{ + BL_PROFILE("LBM::generate_voxel_cracks()"); + + amrex::ParmParse pp("voxel_cracks"); + + // Get grid dimensions from domain + const amrex::Box& domain = geom.Domain(); + const int nx = domain.length(0); + const int ny = domain.length(1); + const int nz = domain.length(2); + + amrex::Print() << "Loading voxel cracks for domain: " << nx << " x " << ny + << " x " << nz << std::endl; + + // Get binary crack file path + std::string crack_file; + if (pp.query("crack_file", crack_file) == 0) { + // Default filename pattern matching mainCrackGenerator.C output + crack_file = "microstructure_nX" + std::to_string(nx) + "_nY" + + std::to_string(ny) + "_nZ" + std::to_string(nz) + ".bin"; + } + + // Read crack pattern from file (auto-detect format) + std::vector crack_data = read_crack_file(crack_file, nx, ny, nz); + + // Initialize all cells as SOLID first + is_fluid.setVal(1); + + // Copy crack data to GPU-accessible memory for CUDA builds + amrex::Gpu::DeviceVector d_crack_data(crack_data.size()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, crack_data.begin(), crack_data.end(), + d_crack_data.begin()); + amrex::Gpu::synchronize(); + + // Get raw pointer for device access + auto const* crack_ptr = d_crack_data.data(); + + // Copy crack data to MultiFab using GPU-compatible approach + // Your file stores in k,j,i (z,y,x) order + for (amrex::MFIter mfi(is_fluid); mfi.isValid(); ++mfi) { + const amrex::Box& box = mfi.validbox(); + amrex::Array4 const& is_fluid_arr = is_fluid.array(mfi); + + amrex::ParallelFor( + box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Convert AMReX (i,j,k) to file index (z,y,x order) + int file_index = k * (nx * ny) + j * nx + i; + // Your binary file: 0 = fluid (tubes), 1 = solid + // AMReX m_is_fluid: 0 = solid, 1 = fluid + // So we need to invert the values + is_fluid_arr(i, j, k, 0) = (crack_ptr[file_index] == 0) ? 1 : 0; + }); + } + amrex::Gpu::synchronize(); + + amrex::Print() << "Voxel crack generation complete" << std::endl; +} + } // namespace lbm diff --git a/Source/LBM.H b/Source/LBM.H index 6f05dbb5..1758ee9c 100644 --- a/Source/LBM.H +++ b/Source/LBM.H @@ -304,6 +304,12 @@ private: // Initial temperature. amrex::Real m_initialTemperature = 1.0 / 3.0; + // Body is isothermal. + bool m_bodyIsIsothermal = false; + + // Body temperature. + amrex::Real m_bodyTemperature = 1.0 / 3.0; + // Adiabatic exponent. amrex::Real m_adiabaticExponent = 5.0 / 3.0; diff --git a/Source/LBM.cpp b/Source/LBM.cpp index f6d7d116..3b1a4c9f 100644 --- a/Source/LBM.cpp +++ b/Source/LBM.cpp @@ -60,6 +60,7 @@ LBM::LBM() m_idata_varnames.push_back("is_fluid"); m_idata_varnames.push_back("eb_boundary"); + m_idata_varnames.push_back("eb_fluid_boundary"); for (const auto& vname : m_macrodata_varnames) { m_lbm_varnames.push_back(vname); } @@ -286,6 +287,9 @@ void LBM::read_parameters() pp.query("initial_temperature", m_initialTemperature); + pp.query("body_is_isothermal", m_bodyIsIsothermal); + pp.query("body_temperature", m_bodyTemperature); + pp.query("adiabatic_exponent", m_adiabaticExponent); pp.query("mean_molecular_mass", m_m_bar); @@ -776,6 +780,8 @@ void LBM::relax_f_to_equilibrium(const int lev) amrex::Real nu = m_nu; amrex::Real dt = m_dts[lev]; + const bool body_is_isothermal = m_bodyIsIsothermal; + amrex::ParallelFor( m_f[lev], m_eq[lev].nGrowVect(), constants::N_MICRO_STATES, [=] AMREX_GPU_DEVICE( @@ -799,6 +805,12 @@ void LBM::relax_f_to_equilibrium(const int lev) f_arr(iv, q) += omega * (eq_arr(iv, q) - f_arr(iv, q)); g_arr(iv, q) += omega * (eq_arr_g(iv, q) - g_arr(iv, q)); + + if (body_is_isothermal) { + if (is_fluid_arrs[nbx](iv, 2) == 1) { + g_arr(iv, q) = eq_arr_g(iv, q); + } + } } }); amrex::Gpu::synchronize(); @@ -818,6 +830,9 @@ void LBM::f_to_macrodata(const int lev) amrex::Real specific_gas_constant = m_R_u / m_m_bar; amrex::Real cv = specific_gas_constant / (m_adiabaticExponent - 1.0); + const bool body_is_isothermal = m_bodyIsIsothermal; + const amrex::Real body_temperature = m_bodyTemperature; + const stencil::Stencil stencil; const auto& evs = stencil.evs; amrex::ParallelFor( @@ -888,6 +903,12 @@ void LBM::f_to_macrodata(const int lev) amrex::Real temperature = get_temperature(two_rho_e, rho, u, v, w, cv); + if (body_is_isothermal) { + if (is_fluid_arrs[nbx](iv, 2) == 1) { + temperature = body_temperature; + } + } + md_arr(iv, constants::TEMPERATURE_IDX) = temperature; md_arr(iv, constants::Q_CORR_X_IDX) = @@ -1135,11 +1156,11 @@ void LBM::MakeNewLevelFromCoarse( f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - - compute_q_corrections(lev); } // Make a new level from scratch using provided BoxArray and @@ -1191,11 +1212,11 @@ void LBM::MakeNewLevelFromScratch( f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - - compute_q_corrections(lev); } void LBM::initialize_f(const int lev) @@ -1258,6 +1279,29 @@ void LBM::initialize_is_fluid(const int lev) } }); + // Compute the boundary cells on the fluid side + const stencil::Stencil stencil; + const auto& evs = stencil.evs; + amrex::ParallelFor( + m_is_fluid[lev], m_is_fluid[lev].nGrowVect() - 1, + [=] AMREX_GPU_DEVICE( + int nbx, int i, int j, int AMREX_D_PICK(, /*k*/, k)) noexcept { + const amrex::IntVect iv(AMREX_D_DECL(i, j, k)); + const auto if_arr = is_fluid_arrs[nbx]; + + bool all_covered = true; + for (int idir = 0; idir < constants::N_MICRO_STATES; idir++) { + const auto& dimvec = evs[idir]; + all_covered &= (if_arr(iv - dimvec, 0) == 1); + } + + if ((all_covered) || (if_arr(iv, 0) == 0)) { + if_arr(iv, 2) = 0; + } else { + if_arr(iv, 2) = 1; + } + }); + m_is_fluid[lev].FillBoundary(Geom(lev).periodicity()); } @@ -1353,12 +1397,12 @@ void LBM::RemakeLevel( f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - compute_q_corrections(lev); - m_ts_new[lev] = time; m_ts_old[lev] = time - constants::SMALL_NUM; } @@ -1907,11 +1951,11 @@ void LBM::read_checkpoint_file() f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - - compute_q_corrections(lev); } } diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index f0ddd600..f22bfa9c 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -95,7 +95,7 @@ endif() if(MARBLES_DIM EQUAL 3) add_test_r(single_cylinder) add_test_r(channel_cylinder) - # add_test_r(channel_cylinder_amr) + add_test_r(cylinder_turek_2d2) endif() add_test_re(sod) diff --git a/Tests/test_files/cylinder_turek_2d2_amr/cylinder_turek_2d2_amr.inp b/Tests/test_files/cylinder_turek_2d2_amr/cylinder_turek_2d2_amr.inp new file mode 100644 index 00000000..3a444e9a --- /dev/null +++ b/Tests/test_files/cylinder_turek_2d2_amr/cylinder_turek_2d2_amr.inp @@ -0,0 +1,65 @@ +max_step = 100000 + +# geometry parameters +geometry.prob_lo = 0.0 0.0 -2.0 +geometry.prob_hi = 660.0 123.0 2.0 +geometry.is_periodic = 0 0 1 + +# timestepping +amr.n_cell = 660 123 4 +amr.max_level = 1 +amr.max_grid_size = 32 +amr.blocking_factor_x = 2 +amr.blocking_factor_y = 1 +amr.blocking_factor_z = 1 +amr.plot_int = 100 +amr.chk_int = 1000 +amr.file_name_digits = 5 + +lbm.bc_lo = 2 1 0 +lbm.bc_hi = 5 1 0 +lbm.dx_outer = 1.0 +lbm.dt_outer = 1.0 + +#nu = u D / Re +#nu = M a D / Re +#nu = M sqrt(gamma R T) D / Re +#nu = 0.1 sqrt(1.4 1 0.03333) 30 / 100 + +lbm.nu = 0.00648 +lbm.alpha = 0.00648 +lbm.initial_temperature = 0.03333 +lbm.adiabatic_exponent = 1.4 +lbm.save_streaming = 0 +lbm.compute_forces = 1 + +lbm.velocity_bc_type = "parabolic" +velocity_bc_parabolic.Mach_ref = 0.1 +velocity_bc_parabolic.normal_dir = 1 +velocity_bc_parabolic.tangential_dir = 0 +velocity_bc_parabolic.initial_density = 1.0 +velocity_bc_parabolic.initial_temperature = 0.03333 +velocity_bc_parabolic.adiabatic_exponent = 1.4 + +lbm.ic_type = "constant" +ic_constant.mach_components = 0.1 0.0 0.0 +ic_constant.density = 1.0 +ic_constant.initial_temperature = 0.03333 +ic_constant.adiabatic_exponent = 1.4 + +eb2.geom_type = "cylinder" +eb2.cylinder_radius = 15.0 +eb2.cylinder_center = 60.0 60.0 0.0 +eb2.cylinder_has_fluid_inside = 0 +eb2.cylinder_height = 1000.0 +eb2.cylinder_direction = 2 + +tagging.refinement_indicators = box +tagging.box.in_box_lo = 30.0 30.0 -2.0 +tagging.box.in_box_hi = 90.0 90.0 2.0 + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 +amrex.the_arena_is_managed = 0 +amrex.abort_on_out_of_gpu_memory = 1 \ No newline at end of file diff --git a/Tests/test_files/isothermal_cracks/isothermal_cracks.inp b/Tests/test_files/isothermal_cracks/isothermal_cracks.inp new file mode 100644 index 00000000..071c7d5d --- /dev/null +++ b/Tests/test_files/isothermal_cracks/isothermal_cracks.inp @@ -0,0 +1,78 @@ +max_step = 1000 + +# geometry parameters - coarse domain +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 60.0 40.0 30.0 +geometry.is_periodic = 0 1 1 + +# timestepping - coarse grid for thick cracks +amr.n_cell = 60 40 30 +amr.max_level = 0 +amr.max_grid_size = 64 +amr.blocking_factor_x = 2 +amr.blocking_factor_y = 2 +amr.blocking_factor_z = 2 +amr.plot_int = 10 +amr.chk_int = 1000 +amr.file_name_digits = 5 + +lbm.bc_lo = 2 0 0 +lbm.bc_hi = 5 0 0 +lbm.dx_outer = 1.0 +lbm.dt_outer = 1.0 + +#Parameters and unit conversion +#For air flowing through 6 millimeter crack at 0.1 meters/second +#Re = D u / nu = 6e-3 * 0.1 / 1.5e-5 = 40 + +#Calculating nu in LB units +# dx = 1 = 1 millimeter. D = 6. Verify in mainCrackGenerator.C that std::vector diameterVector{7, 6, 5} +#nu = D u / Re +#nu = D M a / Re +#nu = D M sqrt(gamma R T) / Re +#Since we are using a compressible code for incompressible problem, Mach can be used as a free parameter +#Setting M=0.1 +#Setting T = 1/30 . Now corresponds to physical temperature at viscosity 1.5e-5 m^2/s, i.e. T=288 K +#Except in combustion, only ratios of temeperature matter. E.g., T = 2/30 now corresponds to 576 K +#R has be set in the code as a ratio of R_air, i.e. R=1 corresponds to 287 J/(kg K) +#Adiabatic exponent for air, gamma=1.4 +#nu = 6 0.1 sqrt(1.4 1 0.03333) / 40 +#alpha = nu / Prandtl. Prandtl number is approx 0.7 for air + +lbm.nu = 3.238e-3 +lbm.alpha = 4.626e-3 +lbm.initial_temperature = 0.03333 +lbm.adiabatic_exponent = 1.4 +lbm.save_streaming = 0 +# lbm.compute_forces = 1 + +lbm.body_is_isothermal = 1 +lbm.body_temperature = 0.066 + +lbm.velocity_bc_type = "parabolic" +velocity_bc_parabolic.Mach_ref = 0.1 +velocity_bc_parabolic.normal_dir = 1 +velocity_bc_parabolic.tangential_dir = 0 +velocity_bc_parabolic.initial_density = 1.0 +velocity_bc_parabolic.initial_temperature = 0.03333 +velocity_bc_parabolic.adiabatic_exponent = 1.4 + +lbm.ic_type = "constant" +ic_constant.mach_components = 0.0 0.0 0.0 +ic_constant.density = 1.0 +ic_constant.initial_temperature = 0.03333 +ic_constant.adiabatic_exponent = 1.4 + +# Use all_regular geometry with voxel crack generation +eb2.geom_type = "all_regular" +eb2.use_voxel_cracks = 1 + +# Voxel crack generation from binary file (also supports .csv) +# Geometry generation codes and tools are in ../../../Tools/isothermal_cracks +voxel_cracks.crack_file = "microstructure_nX60_nY40_nZ30.bin" + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 +amrex.the_arena_is_managed = 0 +amrex.abort_on_out_of_gpu_memory = 1 \ No newline at end of file diff --git a/Tests/test_files/isothermal_cracks/microstructure_nX60_nY40_nZ30.bin b/Tests/test_files/isothermal_cracks/microstructure_nX60_nY40_nZ30.bin new file mode 100644 index 00000000..efcd6348 Binary files /dev/null and b/Tests/test_files/isothermal_cracks/microstructure_nX60_nY40_nZ30.bin differ diff --git a/Tests/test_files/isothermal_cylinder/isothermal_cylinder.inp b/Tests/test_files/isothermal_cylinder/isothermal_cylinder.inp new file mode 100644 index 00000000..5ae5169b --- /dev/null +++ b/Tests/test_files/isothermal_cylinder/isothermal_cylinder.inp @@ -0,0 +1,64 @@ +max_step = 100000 + +# geometry parameters +geometry.prob_lo = 0.0 0.0 -2.0 +geometry.prob_hi = 660.0 123.0 2.0 +geometry.is_periodic = 0 0 1 + +# timestepping +amr.n_cell = 660 123 4 +amr.max_level = 0 +amr.max_grid_size = 32 +amr.blocking_factor_x = 2 +amr.blocking_factor_y = 1 +amr.blocking_factor_z = 1 +amr.plot_int = 100 +amr.chk_int = 1000 +amr.file_name_digits = 5 + +lbm.bc_lo = 2 1 0 +lbm.bc_hi = 5 1 0 +lbm.dx_outer = 1.0 +lbm.dt_outer = 1.0 + +#nu = u D / Re +#nu = M a D / Re +#nu = M sqrt(gamma R T) D / Re +#nu = 0.1 sqrt(1.4 1 0.03333) 30 / 100 + +lbm.nu = 0.00648 +lbm.alpha = 0.00648 +lbm.initial_temperature = 0.03333 +lbm.adiabatic_exponent = 1.4 +lbm.save_streaming = 0 +lbm.compute_forces = 1 + +lbm.body_is_isothermal = 1 +lbm.body_temperature = 0.060 + +lbm.velocity_bc_type = "parabolic" +velocity_bc_parabolic.Mach_ref = 0.1 +velocity_bc_parabolic.normal_dir = 1 +velocity_bc_parabolic.tangential_dir = 0 +velocity_bc_parabolic.initial_density = 1.0 +velocity_bc_parabolic.initial_temperature = 0.03333 +velocity_bc_parabolic.adiabatic_exponent = 1.4 + +lbm.ic_type = "constant" +ic_constant.mach_components = 0.1 0.0 0.0 +ic_constant.density = 1.0 +ic_constant.initial_temperature = 0.03333 +ic_constant.adiabatic_exponent = 1.4 + +eb2.geom_type = "cylinder" +eb2.cylinder_radius = 15.0 +eb2.cylinder_center = 60.0 60.0 0.0 +eb2.cylinder_has_fluid_inside = 0 +eb2.cylinder_height = 1000.0 +eb2.cylinder_direction = 2 + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 +amrex.the_arena_is_managed = 0 +amrex.abort_on_out_of_gpu_memory = 1 \ No newline at end of file diff --git a/Tests/test_files/isothermal_plates/isothermal_plates.inp b/Tests/test_files/isothermal_plates/isothermal_plates.inp new file mode 100644 index 00000000..696bf4b3 --- /dev/null +++ b/Tests/test_files/isothermal_plates/isothermal_plates.inp @@ -0,0 +1,68 @@ +max_step = 100000 + +# geometry parameters +geometry.prob_lo = 0.0 0.0 -2.0 +geometry.prob_hi = 660.0 123.0 2.0 +geometry.is_periodic = 0 0 1 + +# timestepping +amr.n_cell = 660 123 4 +amr.max_level = 0 +amr.max_grid_size = 32 +amr.blocking_factor_x = 2 +amr.blocking_factor_y = 1 +amr.blocking_factor_z = 1 +amr.plot_int = 100 +amr.chk_int = 1000 +amr.file_name_digits = 5 + +lbm.bc_lo = 2 1 0 +lbm.bc_hi = 5 1 0 +lbm.dx_outer = 1.0 +lbm.dt_outer = 1.0 + +#nu = u D / Re +#nu = M a D / Re +#nu = M sqrt(gamma R T) D / Re +#nu = 0.1 sqrt(1.4 1 0.03333) 30 / 100 + +lbm.nu = 0.00648 +lbm.alpha = 0.00648 +lbm.initial_temperature = 0.03333 +lbm.adiabatic_exponent = 1.4 +lbm.save_streaming = 0 +# lbm.compute_forces = 1 + +lbm.body_is_isothermal = 1 +lbm.body_temperature = 0.060 + +lbm.velocity_bc_type = "parabolic" +velocity_bc_parabolic.Mach_ref = 0.1 +velocity_bc_parabolic.normal_dir = 1 +velocity_bc_parabolic.tangential_dir = 0 +velocity_bc_parabolic.initial_density = 1.0 +velocity_bc_parabolic.initial_temperature = 0.03333 +velocity_bc_parabolic.adiabatic_exponent = 1.4 + +lbm.ic_type = "constant" +ic_constant.mach_components = 0.1 0.0 0.0 +ic_constant.density = 1.0 +ic_constant.initial_temperature = 0.03333 +ic_constant.adiabatic_exponent = 1.4 + +# eb2.geom_type = "cylinder" +# eb2.cylinder_radius = 15.0 +# eb2.cylinder_center = 60.0 60.0 0.0 +# eb2.cylinder_has_fluid_inside = 0 +# eb2.cylinder_height = 1000.0 +# eb2.cylinder_direction = 2 + +eb2.geom_type = "plane" +eb2.plane_point = 0 1 0 +eb2.plane_normal = 0.0 -1.0 0 + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 +amrex.the_arena_is_managed = 0 +amrex.abort_on_out_of_gpu_memory = 1 \ No newline at end of file diff --git a/Tools/isothermal_cracks/binaryToTiff.py b/Tools/isothermal_cracks/binaryToTiff.py new file mode 100644 index 00000000..8b4bcb1c --- /dev/null +++ b/Tools/isothermal_cracks/binaryToTiff.py @@ -0,0 +1,113 @@ +import sys +import time +from pathlib import Path + +def import_required_modules(): + """ + Import required dependencies and fail fast with actionable guidance if they + are missing. The caller can install them explicitly (e.g. via `uv pip` + or their preferred environment manager). + """ + + try: + import numpy as np # type: ignore + except ImportError as exc: # pragma: no cover - defensive + raise RuntimeError( + "Missing dependency 'numpy'. Install it via your environment manager," + " e.g. `uv pip install numpy`." + ) from exc + + try: + import tifffile as tf # type: ignore + except ImportError as exc: # pragma: no cover - defensive + raise RuntimeError( + "Missing dependency 'tifffile'. Install it via your environment manager," + " e.g. `uv pip install tifffile`." + ) from exc + + try: + import imagecodecs # type: ignore + compression_available = True + except ImportError: + compression_available = False + + return np, tf, compression_available + + +def load_dimensions(dimensions_path: str): + """Read the domain dimensions from the helper file.""" + + with open(dimensions_path, "r", encoding="utf-8") as file_obj: + dims = file_obj.read().strip().split() + if len(dims) != 3: + raise ValueError( + f"Unexpected contents in {dimensions_path!r}: expected three integers" + ) + return tuple(int(value) for value in dims) + + +def find_binary_file(n_x: int, n_y: int, n_z: int) -> Path: + """Locate the binary file, allowing for optional suffixes like _seedX.""" + + results_dir = Path("results") + pattern = f"microstructure_nX{n_x}_nY{n_y}_nZ{n_z}*.bin" + matches = sorted(results_dir.glob(pattern)) + + if not matches: + raise FileNotFoundError( + f"No binary file matching pattern {pattern!r} in {results_dir.resolve()}" + ) + + if len(matches) > 1: + chosen = matches[0] + print( + "Found multiple matching binary files, using the first one:" + f" {chosen.name}\nCandidates: {[path.name for path in matches]}" + ) + return chosen + + return matches[0] + + +def main() -> int: + try: + np, tf, compression_available = import_required_modules() + except RuntimeError as exc: # pragma: no cover - defensive + print(exc, file=sys.stderr) + return 1 + + print("Starting binary to TIFF conversion...") + start_time = time.time() + + n_x, n_y, n_z = load_dimensions("results/dimensions.txt") + print(f"Dimensions: {n_x} x {n_y} x {n_z}") + + print("Loading binary data...") + binary_path = find_binary_file(n_x, n_y, n_z) + image_array = np.fromfile(str(binary_path), dtype=np.uint16) + image_array = image_array.reshape((n_z, n_y, n_x)) + + print("Saving TIFF...") + tiff_path = binary_path.with_suffix(".tif") + compression_args = {"compression": "lzw"} if compression_available else {} + + try: + tf.imwrite(str(tiff_path), image_array, **compression_args) + if compression_available: + print("Saved with LZW compression") + else: + print("Saved without compression") + except Exception as exc: # pragma: no cover - defensive + print(f"Compression failed: {exc}") + print("Retrying without compression...") + tf.imwrite(str(tiff_path), image_array) + print("Saved without compression") + + end_time = time.time() + print(f"TIFF conversion completed in {end_time - start_time:.2f} seconds") + print(f"TIFF file saved: {tiff_path}") + return 0 + + +if __name__ == "__main__": # pragma: no cover - script entry point + sys.exit(main()) \ No newline at end of file diff --git a/Tools/isothermal_cracks/mainCrackGenerator.C b/Tools/isothermal_cracks/mainCrackGenerator.C new file mode 100644 index 00000000..4227ce15 --- /dev/null +++ b/Tools/isothermal_cracks/mainCrackGenerator.C @@ -0,0 +1,476 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float myReal; + +int main(int argc, char* argv[]) +{ + // Get geometry seed from command line argument or use default + int geometrySeed = 12345; // default seed + if (argc > 1) { + geometrySeed = std::atoi(argv[1]); + } + + std::cout << "Using geometry seed: " << geometrySeed << std::endl; + + const int nX(60), nY(40), nZ(30); // grid + + std::vector diameterVector{7, 6, 5}; // diameter of tubes + + std::vector volumeFractions{ + 0.2, 0.3, 0.4}; // relativevolume fractions of corresponding tubes + // Need not sum up to 1, will be automatically rescaled to 1 + + myReal clusteringFactor = + 0.001; // reduce for more clustering, increase for less clustering + + myReal tubeVolumeFraction = + 0.15; // total volume fraction of all tubes, i.e. the empty space + + bool forMARBLES_LBM_code = + true; // makes tubes 0 and rest of the domain of solid space is 1. + // Otherwise, tubes are >0 and solid space is 0 + + // Scale the volume fractions to a total of 1.0 + myReal sumVF = 0.0; + for (const auto& vf : volumeFractions) sumVF += vf; + for (auto& vf : volumeFractions) + vf = vf / sumVF; // rescale so that the sum is 1. + + for (size_t i = 0; i < volumeFractions.size(); ++i) + volumeFractions[i] *= tubeVolumeFraction; + + // no user input below this line + std::vector radiusVector(diameterVector.size()); + for (size_t i = 0; i < diameterVector.size(); ++i) + radiusVector[i] = diameterVector[i] / 2.0; // convert diameter to radius + + const long maxNumberOfSpheres = static_cast(nX) * nY * nZ; + int maxRadius = + std::ceil(*std::max_element(radiusVector.begin(), radiusVector.end())); + const int pX(maxRadius), pY(maxRadius), pZ(maxRadius); // dummy points + const int nsX(nX + 2 * pX), nsY(nY + 2 * pY), nsZ(nZ + 2 * pZ); + + // Use a more memory-efficient 3D structure instead of 4D vector + std::vector>> grid( + nsZ, std::vector>(nsY, std::vector(nsX, 0))); + + std::default_random_engine generator(geometrySeed); + std::uniform_int_distribution distributionX(pX, pX + nX - 1); + std::uniform_int_distribution distributionY(pY, pY + nY - 1); + std::uniform_int_distribution distributionZ(pZ, pZ + nZ - 1); + + std::uniform_int_distribution distributionPathFluctuation(-1, 1); + + std::vector tags(radiusVector.size()); + for (int i = 0; i < radiusVector.size(); i++) + tags[i] = i + 1; // radiusVector[i]; //i+1; + + std::vector count(radiusVector.size() + 1, 0); + std::vector countTarget(radiusVector.size() + 1, 0); + count[0] = static_cast(nX) * nY * nZ; + + myReal porosity = + 1.0 - tubeVolumeFraction; // porosity is 1 - tube volume fraction + + countTarget[0] = porosity * count[0]; + for (int i = 0; i < countTarget.size(); i++) + countTarget[i + 1] = + volumeFractions[i] * static_cast(nX) * nY * nZ; + + // Free pass + for (int entity = 0; entity < std::max(1, int(radiusVector.size() / 1)); + entity++) { + int y, z; // define y and z outside the loop to make tubes along x + // direction + + long maxIterations = static_cast( + static_cast(clusteringFactor) * + static_cast(maxNumberOfSpheres)); + + for (long i = 0; i < maxIterations; i++) { + if (count[entity + 1] < countTarget[entity + 1]) { + + int radius = std::ceil(radiusVector[entity]); + + int x = pX + (i % nX); + + if (x == pX) // new y and z for each sweep along x + { + y = distributionY(generator); + z = distributionZ(generator); + } else { + y += distributionPathFluctuation(generator) * (i % 2); + z += distributionPathFluctuation(generator) * ((i + 1) % 2); + } + + // make y and z periodic with modulo in case they go out of + // bounds + y = (y - pY + nY) % nY + pY; + z = (z - pZ + nZ) % nZ + pZ; + + // check if y and z are within bounds + if ((y < pY) or (y >= pY + nY) or (z < pZ) or (z >= pZ + nZ)) { + + } else { + + if (grid[z][y][x] == 0) { + count[0] -= 1; + grid[z][y][x] = tags[entity]; + count[1 + entity] += 1; + + for (int zz = z - radius; zz <= z + radius; zz++) + for (int yy = y - radius; yy <= y + radius; yy++) + for (int xx = x - radius; xx <= x + radius; + xx++) { + if (grid[zz][yy][xx] == 0) { + if (std::sqrt( + (xx - x) * (xx - x) + + (yy - y) * (yy - y) + + (zz - z) * (zz - z)) <= + radiusVector[entity]) { + count[0] -= 1; + grid[zz][yy][xx] = tags[entity]; + count[1 + entity] += 1; + } + } + } + } + } + } + } + if (count[0] <= countTarget[0]) { + std::cout << "Target porosity reached during structure generation. " + "Stopping and continuing to write to file" + << std::endl; + break; + } + } + + // Pass with proximity check + long maxIterations = static_cast( + static_cast(1.0) * static_cast(maxNumberOfSpheres)); + + for (int entity = 0; entity < radiusVector.size(); entity++) { + int y, z; // define y and z outside the loop to make tubes along x + // direction + + for (long i = 0; i < maxIterations; i++) { + + if (count[entity + 1] < countTarget[entity + 1]) { + + int radius = std::ceil(radiusVector[entity]); + + // x sweeps along the X direction to make multiple tubes in + // passes Algorithm: modulo if i/nX + int x = pX + (i % nX); + + if (i % nX == 0) // new y and z for each sweep along x + { + y = distributionY(generator); + z = distributionZ(generator); + } else { + y += distributionPathFluctuation(generator) * (i % 2); + z += distributionPathFluctuation(generator) * ((i + 1) % 2); + } + + // make y and z periodic with modulo in case they go out of + // bounds + y = (y - pY + nY) % nY + pY; + z = (z - pZ + nZ) % nZ + pZ; + + // check if y and z are within bounds + if ((y < pY) or (y >= pY + nY) or (z < pZ) or (z >= pZ + nZ)) { + + } else { + + if (grid[z][y][x] == 0) { + int connectivityChecked = + 0; // to check if connectivity is checked for this + // point + + for (int zz = z - radius; zz <= z + radius; zz++) + for (int yy = y - radius; yy <= y + radius; yy++) + for (int xx = x - radius; xx <= x + radius; + xx++) { + + if (grid[zz][yy][xx] == 0) { + + if ((grid[z - (zz - z)][y - (yy - y)] + [x - (xx - x)] != 0) or + (connectivityChecked == + 1)) // check if the point is + // connected to any other point + { + + connectivityChecked = 1; + + if (std::sqrt( + (xx - x) * (xx - x) + + (yy - y) * (yy - y) + + (zz - z) * (zz - z)) <= + radiusVector[entity]) { + count[0] -= 1; + grid[zz][yy][xx] = tags[entity]; + count[1 + entity] += 1; + + if (grid[z][y][x] == 0) { + count[0] -= 1; + grid[z][y][x] = + tags[entity]; + count[1 + entity] += 1; + } + } + } + } + } + } + } + } + } + if (count[0] <= countTarget[0]) { + std::cout << "Target porosity reached during structure generation. " + "Stopping and continuing to write to file" + << std::endl; + break; + } + } + + if (count[0] > countTarget[0]) { + std::cout << "Target porosity not reached during structure generation. " + "Continuing to bridge spheres" + << std::endl; + + // Do bridging if target porosity is not reached + for (int entity = radiusVector.size() - 1; entity > 0; entity--) { + if (count[0] <= countTarget[0]) break; + std::cout << "Bridging entity " << entity << " with radius " + << radiusVector[entity] << std::endl; + + for (int i = 0; i < int(0.01 * maxNumberOfSpheres); i++) { + + if (count[0] <= countTarget[0]) { + break; + } + + int radius = std::ceil(radiusVector[entity]); + + int x = distributionX(generator); + int y = distributionY(generator); + int z = distributionZ(generator); + + if (grid[z][y][x] == 0) { + int connectivityChecked = + 0; // to check if connectivity is checked for this point + + for (int zz = z - radius; zz <= z + radius; zz++) + for (int yy = y - radius; yy <= y + radius; yy++) + for (int xx = x - radius; xx <= x + radius; xx++) { + if (grid[zz][yy][xx] == 0) { + if ((grid[z - (zz - z)][y - (yy - y)] + [x - (xx - x)] != 0) or + (connectivityChecked == + 1)) // check if the point is connected + // to any other point + { + connectivityChecked = 1; + + if (std::sqrt( + (xx - x) * (xx - x) + + (yy - y) * (yy - y) + + (zz - z) * (zz - z)) <= + radiusVector[entity]) { + count[0] -= 1; + grid[zz][yy][xx] = tags[entity]; + count[1 + entity] += 1; + + if (grid[z][y][x] == 0) { + count[0] -= 1; + grid[z][y][x] = tags[entity]; + count[1 + entity] += 1; + } + } + } + } + } + } + } + } + } + + if (count[0] > countTarget[0]) { + std::cout << "WARNING!!! Target porosity could not be achieved" + << std::endl; + } else { + std::cout << "Target porosity reached during bridging. Stopping and " + "continuing to write to file" + << std::endl; + } + + // marbles compatibility: make tubes 0 and solid space 1 + if (forMARBLES_LBM_code) { + for (int k = 0; k < grid.size(); k++) + for (int j = 0; j < grid[0].size(); j++) + for (int i = 0; i < grid[0][0].size(); i++) { + // Make tag > 0 as 0 and 0 as 1 + if (grid[k][j][i] > 0) + grid[k][j][i] = 0; + else + grid[k][j][i] = 1; + } + + // Mark the first physical X plane as fluid + for (int k = 0; k < grid.size(); k++) + for (int j = 0; j < grid[0].size(); j++) grid[k][j][pX] = 0; + + // Mark the last physical X plane as fluid + for (int k = 0; k < grid.size(); k++) + for (int j = 0; j < grid[0].size(); j++) + grid[k][j][grid[0][0].size() - 1 - pX] = 0; + } + + std::filesystem::create_directories("results"); + + std::ofstream ofile; + std::string filename = "./results/microstructure_nX" + std::to_string(nX) + + "_nY" + std::to_string(nY) + "_nZ" + + std::to_string(nZ); + + filename += "_seed" + std::to_string(geometrySeed) + ".csv"; + + ofile.open(filename); + + // Build the entire output as a string first + std::string output; + output.reserve(1000000000); // Reserve space for ~1GB string + + output += "X,Y,Z,tag\n"; + + for (int k = pZ; k < grid.size() - pZ; k++) { + std::cout << 100.0 * (k - pZ) / (grid.size() - 1) << " %" << std::endl; + for (int j = pY; j < grid[0].size() - pY; j++) + for (int i = pX; i < grid[0][0].size() - pX; i++) { + output += std::to_string(i - pX) + "," + + std::to_string(j - pY) + "," + + std::to_string(k - pZ) + "," + + std::to_string(grid[k][j][i]) + "\n"; + } + } + + // Write everything at once + std::cout << "Writing to file: " << filename << std::endl; + ofile << output; + ofile.close(); + std::cout << "File written: " << filename << std::endl; + + // Also write binary data for faster TIFF conversion + std::string binFilename = "./results/microstructure_nX" + + std::to_string(nX) + "_nY" + std::to_string(nY) + + "_nZ" + std::to_string(nZ) + "_seed" + + std::to_string(geometrySeed) + ".bin"; + std::ofstream binFile(binFilename, std::ios::binary); + + std::cout << "Writing binary file for TIFF conversion..." << std::endl; + for (int k = pZ; k < grid.size() - pZ; k++) { + if (k % 100 == 0) + std::cout << "Binary write progress: " + << 100.0 * (k - pZ) / (grid.size() - 2 * pZ) << " %" + << std::endl; + for (int j = pY; j < grid[0].size() - pY; j++) + for (int i = pX; i < grid[0][0].size() - pX; i++) { + uint16_t value = static_cast(grid[k][j][i]); + binFile.write( + reinterpret_cast(&value), sizeof(uint16_t)); + } + } + binFile.close(); + std::cout << "Binary file written: " << binFilename << std::endl; + + // Write dimensions file for Python script + std::ofstream dimFile("./results/dimensions.txt"); + dimFile << nX << " " << nY << " " << nZ << std::endl; + dimFile.close(); + + std::cout << "Open with Paraview with CSV Reader > Apply" << std::endl; + std::cout << "option + space > Table to Structured Grid" << std::endl; + std::cout << "Whole Extent: 0 to " << nX - 1 << std::endl; + std::cout << " 0 to " << nY - 1 << std::endl; + std::cout << " 0 to " << nZ - 1 << std::endl; + std::cout << "X Column X" << std::endl; + std::cout << "Y Column Y" << std::endl; + std::cout << "Z Column Z" << std::endl; + + auto run_with_uvx = [](const std::string& script, + const std::string& uvx_packages) -> int { + std::string uvx_command = "uvx " + uvx_packages + " python " + script; + int status = std::system(uvx_command.c_str()); + if (status != 0) { + std::cout << "uvx command failed (status " << status + << "). Falling back to `python " << script << "`." + << std::endl; + status = std::system(("python " + script).c_str()); + if (status != 0) { + std::cout << "`python` executable not available or script " + "failed (status " + << status << "). Attempting `python3` fallback." + << std::endl; + status = std::system(("python3 " + script).c_str()); + } + } + return status; + }; + + const std::string binary_uvx_packages = + "--with numpy --with tifffile --with imagecodecs"; + const std::string cad_uvx_packages = + "--with numpy --with tifffile --with imagecodecs --with scikit-image " + "--with trimesh"; + + if (std::filesystem::exists("binaryToTiff.py")) { + std::cout << "Running Python TIFF conversion (via uvx)..." << std::endl; + int result = run_with_uvx("binaryToTiff.py", binary_uvx_packages); + if (result == 0) { + std::cout << "TIFF conversion completed successfully!" << std::endl; + + if (std::filesystem::exists("tiffToCAD.py")) { + std::cout << "Running Python CAD conversion (via uvx)..." + << std::endl; + int cadResult = run_with_uvx("tiffToCAD.py", cad_uvx_packages); + if (cadResult == 0) { + std::cout << "CAD conversion completed successfully!" + << std::endl; + std::cout << "STL files saved in cad_exports/ directory" + << std::endl; + } else { + std::cout << "CAD conversion failed. Run manually: uvx " + << cad_uvx_packages << " python tiffToCAD.py" + << std::endl; + } + } else { + std::cout << "tiffToCAD.py not found. To convert to CAD " + "formats, run: uvx " + << cad_uvx_packages << " python tiffToCAD.py" + << std::endl; + } + + } else { + std::cout << "TIFF conversion failed. Run manually: uvx " + << binary_uvx_packages << " python binaryToTiff.py" + << std::endl; + } + } else { + std::cout << "Python script not found. To convert to TIFF, run: uvx " + << binary_uvx_packages << " python binaryToTiff.py" + << std::endl; + } + + return 0; +} diff --git a/Tools/isothermal_cracks/tiffToCAD.py b/Tools/isothermal_cracks/tiffToCAD.py new file mode 100644 index 00000000..7fc2c760 --- /dev/null +++ b/Tools/isothermal_cracks/tiffToCAD.py @@ -0,0 +1,401 @@ +import sys +import time +import os +import traceback + + +def import_required_modules(): + """ + Import the Python dependencies required for TIFF to CAD conversion. + + The function raises a RuntimeError with installation guidance if a + dependency is missing so the caller can install it explicitly (for + example via `uv pip install `). + """ + + try: + import numpy as np # type: ignore + except ImportError as exc: # pragma: no cover - defensive + raise RuntimeError( + "Missing dependency 'numpy'. Install it, e.g. `uv pip install numpy`." + ) from exc + + try: + import tifffile as tf # type: ignore + except ImportError as exc: # pragma: no cover - defensive + raise RuntimeError( + "Missing dependency 'tifffile'. Install it, e.g. `uv pip install tifffile`." + ) from exc + + try: + from skimage import measure # type: ignore + except ImportError as exc: # pragma: no cover - defensive + raise RuntimeError( + "Missing dependency 'scikit-image' (provides skimage.measure)." + " Install it, e.g. `uv pip install scikit-image`." + ) from exc + + try: + import trimesh # type: ignore + except ImportError as exc: # pragma: no cover - defensive + raise RuntimeError( + "Missing dependency 'trimesh'. Install it, e.g. `uv pip install trimesh`." + ) from exc + + try: + import FreeCAD # type: ignore + import Part # type: ignore + import Mesh # type: ignore + freecad_available = True + except ImportError: + FreeCAD = None + Part = None + Mesh = None + freecad_available = False + + return { + "np": np, + "tf": tf, + "measure": measure, + "trimesh": trimesh, + "freecad_available": freecad_available, + "FreeCAD": FreeCAD, + "Part": Part, + "Mesh": Mesh, + } + + +def find_tiff_file(): + """Find the TIFF file in results directory or current directory""" + # Check current directory first + current_dir = "." + tiff_files = [f for f in os.listdir(current_dir) if f.endswith('.tif')] + + if tiff_files: + tiff_file = tiff_files[0] + if len(tiff_files) > 1: + print(f"Multiple TIFF files found in current directory: {tiff_files}") + print(f"Using: {tiff_file}") + return tiff_file + + # Check results directory + results_dir = "results" + if os.path.exists(results_dir): + tiff_files = [f for f in os.listdir(results_dir) if f.endswith('.tif')] + if tiff_files: + tiff_file = os.path.join(results_dir, tiff_files[0]) + if len(tiff_files) > 1: + print(f"Multiple TIFF files found in results directory: {tiff_files}") + print(f"Using: {tiff_files[0]}") + return tiff_file + + raise FileNotFoundError("No TIFF files found in current directory or results/") + +def generate_mesh_for_tag( + volume_data, + tag_value, + np_module, + measure_module, + trimesh_module, + smoothing=True, +): + """Generate mesh for a specific tag value""" + if tag_value == -1: + # Special case for combined mesh + binary_mask = (volume_data > 0).astype(np_module.uint8) + print(f"Generating combined mesh for all solid materials...") + else: + print(f"Generating mesh for tag {tag_value}...") + binary_mask = (volume_data == tag_value).astype(np_module.uint8) + + # Check if tag exists + if np_module.sum(binary_mask) == 0: + print(f"Warning: No voxels found for tag {tag_value if tag_value != -1 else 'combined'}") + return None + + print( + f"Found {np_module.sum(binary_mask)} voxels for tag " + f"{tag_value if tag_value != -1 else 'combined'}" + ) + + # Generate mesh using marching cubes + try: + verts, faces, normals, values = measure_module.marching_cubes( + binary_mask, + level=0.5, + spacing=(1.0, 1.0, 1.0) + ) + + print(f"Generated mesh: {len(verts)} vertices, {len(faces)} faces") + + # Create trimesh object + mesh = trimesh_module.Trimesh(vertices=verts, faces=faces, vertex_normals=normals) + + # Apply smoothing using correct trimesh API + if smoothing and len(verts) > 0: + try: + # Use smoothing filter if available + mesh = mesh.smoothed() + except AttributeError: + try: + # Alternative smoothing method + mesh = trimesh_module.smoothing.filter_laplacian(mesh, iterations=1) + except (AttributeError, ImportError): + # Skip smoothing if not available + print("Smoothing not available - using original mesh") + + # Clean up mesh + mesh.remove_degenerate_faces() + mesh.remove_duplicate_faces() + mesh.remove_unreferenced_vertices() + + if not mesh.is_empty: + print(f"Final mesh: {len(mesh.vertices)} vertices, {len(mesh.faces)} faces") + return mesh + else: + print(f"Warning: Empty mesh generated for tag {tag_value}") + return None + + except Exception as e: + print(f"Error generating mesh for tag {tag_value}: {e}") + return None + +def export_stl(mesh, filename, tag_value): + """Export mesh to STL format""" + try: + if tag_value == -1: + stl_filename = f"{filename}_combined.stl" + else: + stl_filename = f"{filename}_tag{tag_value}.stl" + mesh.export(stl_filename) + print(f"Exported STL: {stl_filename}") + return stl_filename + except Exception as e: + print(f"Error exporting STL for tag {tag_value}: {e}") + return None + +def export_step_iges_freecad( + mesh, + filename, + tag_value, + freecad_modules, + format='step', +): + """Export mesh to STEP/IGES using FreeCAD""" + FreeCAD = freecad_modules[0] + Part = freecad_modules[1] + Mesh_module = freecad_modules[2] + + if FreeCAD is None or Part is None or Mesh_module is None: + print(f"FreeCAD not available - cannot export {format.upper()}") + return None + + try: + # Create FreeCAD document + doc = FreeCAD.newDocument() + + # Convert trimesh to FreeCAD mesh + freecad_mesh = Mesh_module.Mesh() + + # Add triangles to FreeCAD mesh + for face in mesh.faces: + triangle = [mesh.vertices[face[0]], mesh.vertices[face[1]], mesh.vertices[face[2]]] + freecad_mesh.addFacet(triangle[0], triangle[1], triangle[2]) + + # Create mesh object in document + mesh_obj = doc.addObject("Mesh::Feature", f"Tag{tag_value}") + mesh_obj.Mesh = freecad_mesh + + # Export based on format + if format.lower() == 'step': + if tag_value == -1: + output_filename = f"{filename}_combined.step" + else: + output_filename = f"{filename}_tag{tag_value}.step" + Part.export([mesh_obj], output_filename) + + elif format.lower() == 'iges': + if tag_value == -1: + output_filename = f"{filename}_combined.iges" + else: + output_filename = f"{filename}_tag{tag_value}.iges" + Part.export([mesh_obj], output_filename) + else: + raise ValueError(f"Unsupported export format: {format}") + + # Close document + FreeCAD.closeDocument(doc.Name) + + print(f"Exported {format.upper()}: {output_filename}") + return output_filename + + except Exception as e: + print(f"Error exporting {format.upper()} for tag {tag_value}: {e}") + return None + +def main() -> int: + print("Starting TIFF to CAD conversion...") + start_time = time.time() + + try: + modules = import_required_modules() + except RuntimeError as exc: # pragma: no cover - defensive + print(exc, file=sys.stderr) + return 1 + + np_module = modules["np"] + tf_module = modules["tf"] + measure_module = modules["measure"] + trimesh_module = modules["trimesh"] + freecad_available = modules["freecad_available"] + freecad_modules = ( + modules["FreeCAD"], + modules["Part"], + modules["Mesh"], + ) + + if freecad_available: + print("FreeCAD available for STEP/IGES export") + else: + print("FreeCAD not available - only STL export supported") + + try: + tiff_file = find_tiff_file() + print(f"Loading TIFF file: {tiff_file}") + + volume_data = tf_module.imread(tiff_file) + print(f"TIFF dimensions: {volume_data.shape}") + print(f"Data type: {volume_data.dtype}") + + unique_tags = np_module.unique(volume_data) + unique_tags = unique_tags[unique_tags > 0] + print(f"Found tags: {unique_tags}") + + if len(unique_tags) == 0: + print("No non-zero tags found in TIFF file") + return 0 + + output_dir = "cad_exports" + os.makedirs(output_dir, exist_ok=True) + + base_filename = os.path.splitext(os.path.basename(tiff_file))[0] + + exported_files = [] + + max_tags_to_process = min(100, len(unique_tags)) + print( + f"Processing first {max_tags_to_process} tags (out of {len(unique_tags)})" + " for faster testing..." + ) + + for index, tag in enumerate(unique_tags[:max_tags_to_process]): + print(f"\nProcessing tag {tag} ({index + 1}/{max_tags_to_process})...") + + mesh = generate_mesh_for_tag( + volume_data, + tag, + np_module, + measure_module, + trimesh_module, + smoothing=False, + ) + + if mesh is not None: + output_base = os.path.join(output_dir, base_filename) + + stl_file = export_stl(mesh, output_base, tag) + if stl_file: + exported_files.append(stl_file) + + if freecad_available: + step_file = export_step_iges_freecad( + mesh, + output_base, + tag, + freecad_modules, + "step", + ) + if step_file: + exported_files.append(step_file) + + iges_file = export_step_iges_freecad( + mesh, + output_base, + tag, + freecad_modules, + "iges", + ) + if iges_file: + exported_files.append(iges_file) + else: + print("For STEP/IGES export, install FreeCAD (e.g. via conda).") + + print(f"\nGenerating combined mesh...") + combined_mesh = generate_mesh_for_tag( + volume_data, + -1, + np_module, + measure_module, + trimesh_module, + smoothing=False, + ) + + if combined_mesh is not None: + output_base = os.path.join(output_dir, base_filename) + + stl_file = export_stl(combined_mesh, output_base, -1) + if stl_file: + exported_files.append(stl_file) + + if freecad_available: + step_file = export_step_iges_freecad( + combined_mesh, + output_base, + -1, + freecad_modules, + "step", + ) + if step_file: + exported_files.append(step_file) + + iges_file = export_step_iges_freecad( + combined_mesh, + output_base, + -1, + freecad_modules, + "iges", + ) + if iges_file: + exported_files.append(iges_file) + + end_time = time.time() + print(f"\nCAD conversion completed in {end_time - start_time:.2f} seconds") + print(f"Exported {len(exported_files)} files:") + for file_path in exported_files: + if os.path.exists(file_path): + file_size = os.path.getsize(file_path) / (1024 * 1024) + print(f" - {file_path} ({file_size:.1f} MB)") + else: + print(f" - {file_path} (file not found)") + + print(f"\nFiles saved in: {output_dir}/") + print("\nCAD files can be opened in:") + print("- STL: Blender, MeshLab, FreeCAD, SolidWorks, Fusion 360") + print("- STEP: FreeCAD, SolidWorks, Fusion 360, Onshape") + print("- IGES: FreeCAD, SolidWorks, Fusion 360") + + if len(unique_tags) > max_tags_to_process: + print( + f"\nNote: Only processed {max_tags_to_process} out of {len(unique_tags)} tags." + ) + print("To process all tags, change max_tags_to_process in the script.") + + return 0 + + except Exception as exc: # pragma: no cover - defensive + print(f"Error during conversion: {exc}") + traceback.print_exc() + return 1 + +if __name__ == "__main__": # pragma: no cover - script entry point + sys.exit(main()) \ No newline at end of file