# Python script for writing input file for C++ LBM-simulation
# -----------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from skimage import io

# ========================================== 
# Input:
# ------

# LBM computation domain size:
# (Should correspond to the size of the input tif image.)
lx = 200
ly = 200
lz = 200

# 3D binary tif image to be read:
# (Should be in the same folder as this script.)
my_tif_image = 'Sand_002_00_solid_phase_400x400x400_b_bin2.tif'

# ========================================== 
# Functions:
# ----------
 
def read_obstacles():
# Read LBM geometry data from binarised tif image (1 = solid, 0 = air)

    im = io.imread(my_tif_image)
    
    obst = np.zeros((lx,ly,lz)) # matrix including the solid obstacles

    obst = im
                    
    obst = obst.astype(int)
    return obst
#===========================================
def write_obst_file(obst):
# Write the input file    
    my_filename = 'geometry_data.txt'

    solid_count = 0
    pore_count = 0

    with open(my_filename, "a") as log:

        for z in range (0, lz):
            for y in range (0, ly):
                for x in range (0, lx): 
                    log.write("{0}\n".format(str(obst[x,y,z])))
                    
                    if (obst[x,y,z] == 1):
                        solid_count += 1
                    if (obst[x,y,z] == 0):
                        pore_count += 1
                            
    all_cells = pore_count + solid_count
    porosity = pore_count/all_cells
    void_ratio = porosity/(1-porosity)

    return all_cells,pore_count,solid_count,porosity,void_ratio

# ========================================== 
# Main program:
# -------------

print("Reading obstacles from image file...")
obst = read_obstacles()
print("Done.")

print("===========================================")

print("Writing input file for LBM-simulation...")
all_cells,pore_count,solid_count,porosity,void_ratio = write_obst_file(obst)
print("Done.")

print("===========================================")
print("number of solid cells = ", solid_count)
print("number of pore cells = ", pore_count)
print("===========================================")
print("n = ", porosity)
print("e = ", void_ratio)
print("===========================================")