Navier
@author: Alexandre Sac–Morane alexandre.sac-morane@enpc.fr
Generate a synthetic microstructure with addition of the word NAVIER.
A similar consists in using these shape to create the pore space.
Expand source code
#-------------------------------------------------------------------------------
# Librairies
#-------------------------------------------------------------------------------
import skfmm, pickle
import numpy as np
import matplotlib.pyplot as plt
import porespy as ps
from scipy.ndimage import label, binary_dilation
#-------------------------------------------------------------------------------
# Functions
#-------------------------------------------------------------------------------
from numpy_to_vtk import write_vtk_structured_points
from fast_marching_tortuosity import compute_tortuosity_fast_marching
from compute_minkowski import compute_minkowski
#-------------------------------------------------------------------------------
# User
#-------------------------------------------------------------------------------
dim_sample = 150 # -
porosity = 0.5 # -
dim_interface = 4 # -
width_tampon = 4 # -
sample_id = 'special'
#-------------------------------------------------------------------------------
# Tampon
#-------------------------------------------------------------------------------
tampon = np.array([
# N A V I E R
[1,0,0,0,0,0,1, 0,0,1,1,1,0,0, 1,0,0,0,0,0,1, 1,1,1,1,1,1,1, 1,1,1,1,1,1,1, 0,1,1,1,1,0],
[1,1,0,0,0,0,1, 0,1,0,0,0,1,0, 1,0,0,0,0,0,1, 0,0,0,1,0,0,0, 1,0,0,0,0,0,0, 1,0,0,0,0,1],
[1,0,1,0,0,0,1, 1,0,0,0,0,0,1, 0,1,0,0,0,1,0, 0,0,0,1,0,0,0, 1,0,0,0,0,0,0, 1,0,0,0,0,1],
[1,0,0,1,0,0,1, 1,1,1,1,1,1,1, 0,0,1,0,1,0,0, 0,0,0,1,0,0,0, 1,1,1,1,1,0,0, 1,1,1,1,1,0],
[1,0,0,0,1,0,1, 1,0,0,0,0,0,1, 0,0,1,0,1,0,0, 0,0,0,1,0,0,0, 1,0,0,0,0,0,0, 1,0,0,1,0,0],
[1,0,0,0,0,1,1, 1,0,0,0,0,0,1, 0,0,0,1,0,0,0, 0,0,0,1,0,0,0, 1,0,0,0,0,0,0, 1,0,0,0,1,0],
[1,0,0,0,0,0,1, 1,0,0,0,0,0,1, 0,0,0,1,0,0,0, 1,1,1,1,1,1,1, 1,1,1,1,1,1,1, 1,0,0,0,0,1]
])
tampon = np.kron(tampon, np.ones((3, 3), dtype=int))
tampon = binary_dilation(tampon,structure=np.ones((2,2)))
#fig, (ax1) = plt.subplots(1,1, figsize=(16,9),num=1)
#ax1.imshow(tampon_epais, cmap='Greys')
#plt.savefig('todelete.png')
#plt.close()
#-------------------------------------------------------------------------------
# Generate the binary microstructure
#-------------------------------------------------------------------------------
# generate the sample
M_bin = np.zeros((dim_sample, dim_sample, dim_sample))
# init a counter
count = 0
# insert NAVIER randomly until the target porosity is reached
while 1-np.sum(M_bin)/(dim_sample**3) > porosity:
# generate a random cube
position_000 = np.random.rand(3) * dim_sample
count += 1
if count % 500 == 0:
print("Inserted NAVIERs:", count, " Current porosity:", round(1-np.sum(M_bin)/(dim_sample**3),2),'/', porosity)
# create the NAVIER
for i_x in range(int(position_000[0]), int(position_000[0])+tampon.shape[0]):
i_x_tampom = i_x-int(position_000[0])
for i_y in range(int(position_000[1]), int(position_000[1])+tampon.shape[1]):
i_y_tampon = i_y-int(position_000[1])
# determine the value
value = tampon[i_x_tampom, i_y_tampon]
for i_z in range(int(position_000[2]), int(position_000[2])+width_tampon):
# periodic condition (x-axis)
if i_x >= dim_sample :
i_x = i_x - dim_sample
# periodic condition (y-axis)
if i_y >= dim_sample :
i_y = i_y - dim_sample
# periodic condition (z-axis)
if i_z >= dim_sample :
i_z = i_z - dim_sample
# assign
if M_bin[i_x, i_y, i_z] == 0:
M_bin[i_x, i_y, i_z] = value
'''# generate the sample
M_bin = np.zeros((tampon.shape[1]+width_tampon, tampon.shape[1]+width_tampon, tampon.shape[0]))
# faces
for i in range(width_tampon):
M_bin[i, :-width_tampon, :] = np.flip(np.flip(np.transpose(tampon, (1, 0)), 0), 1)
# faces
for i in range(width_tampon):
M_bin[width_tampon:, i, :] = np.flip(np.transpose(tampon, (1, 0)), 1)
# face
for i in range(width_tampon):
M_bin[:-width_tampon, -1-i, :] = np.flip(np.flip(np.transpose(tampon, (1, 0)), 1), 0)
# face
for i in range(width_tampon):
M_bin[-1-i, width_tampon:, :] = np.flip(np.transpose(tampon, (1, 0)), 1)
# blob in the center
M_bin[width_tampon:-width_tampon,width_tampon:-width_tampon,:] = ps.generators.blobs(shape=[M_bin.shape[0]-2*width_tampon, M_bin.shape[1]-2*width_tampon, M_bin.shape[2]], porosity=0.8, blobiness=0.5, periodic=False)
# vtk file
# change the array structure to verify the function
Microstructure_vtk = np.transpose(M_bin, (2, 1, 0))
write_vtk_structured_points('vtk/navier/navier_'+sample_id+'.vtk', Microstructure_vtk, spacing=(1.0, 1.0, 1.0), origin=(0, 0, 0), binary=False)
raise ValueError('stop')'''
#-------------------------------------------------------------------------------
# Compute the sdf
#-------------------------------------------------------------------------------
# Extension of the sample (considering the periodic condition)
# this operation is conducted on the 3 axes
M_bin_extended = np.zeros((dim_sample+2*dim_interface, dim_sample+2*dim_interface, dim_sample+2*dim_interface))
# main sample
M_bin_extended[dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample] = M_bin
# periodicity -x
M_bin_extended[:dim_interface, dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample] = M_bin[-dim_interface:, :, :]
# periodicity +x
M_bin_extended[dim_interface+dim_sample:, dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample] = M_bin[:dim_interface, :, :]
# periodicity -y
M_bin_extended[dim_interface:dim_interface+dim_sample, :dim_interface, dim_interface:dim_interface+dim_sample] = M_bin[:, -dim_interface:, :]
# periodicity +y
M_bin_extended[dim_interface:dim_interface+dim_sample, dim_interface+dim_sample:, dim_interface:dim_interface+dim_sample] = M_bin[:, :dim_interface, :]
# periodicity -z
M_bin_extended[dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample, :dim_interface] = M_bin[:, :, -dim_interface:]
# periodicity +z
M_bin_extended[dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample, dim_interface+dim_sample:] = M_bin[:, :, :dim_interface]
# periodicity -x-y-z
M_bin_extended[:dim_interface, :dim_interface, :dim_interface] = M_bin[-dim_interface:, -dim_interface:, -dim_interface:]
# periodicity -x-y+z
M_bin_extended[:dim_interface, :dim_interface, dim_interface+dim_sample:] = M_bin[-dim_interface:, -dim_interface:, :dim_interface]
# periodicity -x+y-z
M_bin_extended[:dim_interface, dim_interface+dim_sample:, :dim_interface] = M_bin[-dim_interface:, :dim_interface, -dim_interface:]
# periodicity -x+y+z
M_bin_extended[:dim_interface, dim_interface+dim_sample:, dim_interface+dim_sample:] = M_bin[-dim_interface:, :dim_interface, :dim_interface]
# periodicity +x-y-z
M_bin_extended[dim_interface+dim_sample:, :dim_interface, :dim_interface] = M_bin[:dim_interface, -dim_interface:, -dim_interface:]
# periodicity +x-y+z
M_bin_extended[dim_interface+dim_sample:, :dim_interface, dim_interface+dim_sample:] = M_bin[:dim_interface, -dim_interface:, :dim_interface]
# periodicity +x+y-z
M_bin_extended[dim_interface+dim_sample:, dim_interface+dim_sample:, :dim_interface] = M_bin[:dim_interface, :dim_interface, -dim_interface:]
# periodicity +x+y+z
M_bin_extended[dim_interface+dim_sample:, dim_interface+dim_sample:, dim_interface+dim_sample:] = M_bin[:dim_interface, :dim_interface, :dim_interface]
# compute the sdf on the extended sample
M_sd_extended = skfmm.distance(M_bin_extended-0.5, dx = np.array([1, 1, 1]))
# extract the sdf for the original sample
M_sd = M_sd_extended[dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample, dim_interface:dim_interface+dim_sample]
#-------------------------------------------------------------------------------
# Compute the microstructure
#-------------------------------------------------------------------------------
Microstructure = np.zeros((dim_sample, dim_sample, dim_sample))
for i_x in range(dim_sample):
for i_y in range(dim_sample):
for i_z in range(dim_sample):
if M_sd[i_x, i_y, i_z] > dim_interface/2: # inside the grain
Microstructure[i_x, i_y, i_z] = 1
elif M_sd[i_x, i_y, i_z] < -dim_interface/2: # outside the grain
Microstructure[i_x, i_y, i_z] = 0
else : # in the interface
Microstructure[i_x, i_y, i_z] = 0.5 + M_sd[i_x, i_y, i_z]/dim_interface
#-------------------------------------------------------------------------------
# Output
#-------------------------------------------------------------------------------
# plot
#fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2, figsize=(16,9),num=1)
#ax1.imshow(Microstructure[int(dim_sample//2), :, :], cmap='Greys', vmin=0, vmax=1)
#ax1.set_xlabel('Z axis')
#ax1.set_ylabel('Y axis')
#ax2.imshow(Microstructure[:, int(dim_sample//2), :], cmap='Greys', vmin=0, vmax=1)
#ax2.set_xlabel('Z axis')
#ax2.set_ylabel('X axis')
#ax3.imshow(Microstructure[:, :, int(dim_sample//2)], cmap='Greys', vmin=0, vmax=1)
#ax3.set_xlabel('Y axis')
#ax3.set_ylabel('X axis')
#plt.savefig('plot_microstructure.png')
#plt.close()
# save
dict_fft = {'M_microstructure': Microstructure}
with open('fft/navier/dict_fft_'+sample_id, 'wb') as handle:
pickle.dump(dict_fft, handle, protocol=pickle.HIGHEST_PROTOCOL)
# vtk file
# change the array structure to verify the function
Microstructure_vtk = np.transpose(Microstructure, (2, 1, 0))
write_vtk_structured_points('vtk/navier/navier_'+sample_id+'.vtk', Microstructure_vtk, spacing=(1.0, 1.0, 1.0), origin=(0, 0, 0), binary=False)
#-------------------------------------------------------------------------------
# Check the connectivity and extract the connected pores
#-------------------------------------------------------------------------------
# extend the sample (considering the periodic condition)
# this operation is conducted on the 3 axes
M_bin_extended_x = np.zeros((dim_sample, dim_sample+2*dim_sample, dim_sample+2*dim_sample))
for c in range(3):
for h in range(3):
M_bin_extended_x[:, c*dim_sample:(c+1)*dim_sample, h*dim_sample:(h+1)*dim_sample] = M_bin
M_bin_extended_y = np.zeros((dim_sample+2*dim_sample, dim_sample, dim_sample+2*dim_sample))
for l in range(3):
for h in range(3):
M_bin_extended_y[l*dim_sample:(l+1)*dim_sample, :, h*dim_sample:(h+1)*dim_sample] = M_bin
M_bin_extended_z = np.zeros((dim_sample+2*dim_sample, dim_sample+2*dim_sample, dim_sample))
for l in range(3):
for c in range(3):
M_bin_extended_z[l*dim_sample:(l+1)*dim_sample, c*dim_sample:(c+1)*dim_sample, :] = M_bin
# label the pores in the extended samples
labelled_image_x, num_features = label(1-M_bin_extended_x)
labelled_image_y, num_features = label(1-M_bin_extended_y)
labelled_image_z, num_features = label(1-M_bin_extended_z)
# extract sample from the labelled images
extracted_labelled_image_x = labelled_image_x[:, dim_sample:2*dim_sample, dim_sample:2*dim_sample]
extracted_labelled_image_y = labelled_image_y[dim_sample:2*dim_sample, :, dim_sample:2*dim_sample]
extracted_labelled_image_z = labelled_image_z[dim_sample:2*dim_sample, dim_sample:2*dim_sample, :]
# extract the faces
face_xm = extracted_labelled_image_x[0, :, :]
face_xp = extracted_labelled_image_x[-1, :, :]
face_ym = extracted_labelled_image_y[:, 0, :]
face_yp = extracted_labelled_image_y[:, -1, :]
face_zm = extracted_labelled_image_z[:, :, 0]
face_zp = extracted_labelled_image_z[:, :, -1]
# read the unique values on each face (ignore the background 0)
# check if a common id is present on each pair of faces
# create a microstructure with only the connected pores
M_inv_bin_extended_connected_x = np.zeros_like(M_bin_extended_x)
connected_x = False
for id_m in np.unique(face_xm)[1:]:
if id_m in np.unique(face_xp)[1:]:
connected_x = True
M_inv_bin_extended_connected_x = M_inv_bin_extended_connected_x + (labelled_image_x == id_m).astype(int)
M_inv_bin_extended_connected_y = np.zeros_like(M_bin_extended_y)
connected_y = False
for id_m in np.unique(face_ym)[1:]:
if id_m in np.unique(face_yp)[1:]:
connected_y = True
M_inv_bin_extended_connected_y = M_inv_bin_extended_connected_y + (labelled_image_y == id_m).astype(int)
M_inv_bin_extended_connected_z = np.zeros_like(M_bin_extended_z)
connected_z = False
for id_m in np.unique(face_zm)[1:]:
if id_m in np.unique(face_zp)[1:]:
connected_z = True
M_inv_bin_extended_connected_z = M_inv_bin_extended_connected_z + (labelled_image_z == id_m).astype(int)
if not(connected_x and connected_y and connected_z):
print('not connected pores in all directions')
#-------------------------------------------------------------------------------
# Check the connectivity and extract the connected solids
#-------------------------------------------------------------------------------
# label the solids in the extended samples
labelled_image_x, num_features = label(M_bin_extended_x)
labelled_image_y, num_features = label(M_bin_extended_y)
labelled_image_z, num_features = label(M_bin_extended_z)
# extract sample from the labelled images
extracted_labelled_image_x = labelled_image_x[:, dim_sample:2*dim_sample, dim_sample:2*dim_sample]
extracted_labelled_image_y = labelled_image_y[dim_sample:2*dim_sample, :, dim_sample:2*dim_sample]
extracted_labelled_image_z = labelled_image_z[dim_sample:2*dim_sample, dim_sample:2*dim_sample, :]
# extract the faces
face_xm = extracted_labelled_image_x[0, :, :]
face_xp = extracted_labelled_image_x[-1, :, :]
face_ym = extracted_labelled_image_y[:, 0, :]
face_yp = extracted_labelled_image_y[:, -1, :]
face_zm = extracted_labelled_image_z[:, :, 0]
face_zp = extracted_labelled_image_z[:, :, -1]
# read the unique values on each face (ignore the background 0)
# check if a common id is present on each pair of faces
# create a microstructure with only the connected solids
M_bin_extended_connected_x = np.zeros_like(M_bin_extended_x)
connected_x = False
for id_m in np.unique(face_xm)[1:]:
if id_m in np.unique(face_xp)[1:]:
connected_x = True
M_bin_extended_connected_x = M_bin_extended_connected_x + (labelled_image_x == id_m).astype(int)
M_bin_extended_connected_y = np.zeros_like(M_bin_extended_y)
connected_y = False
for id_m in np.unique(face_ym)[1:]:
if id_m in np.unique(face_yp)[1:]:
connected_y = True
M_bin_extended_connected_y = M_bin_extended_connected_y + (labelled_image_y == id_m).astype(int)
M_bin_extended_connected_z = np.zeros_like(M_bin_extended_z)
connected_z = False
for id_m in np.unique(face_zm)[1:]:
if id_m in np.unique(face_zp)[1:]:
connected_z = True
M_bin_extended_connected_z = M_bin_extended_connected_z + (labelled_image_z == id_m).astype(int)
if not(connected_x and connected_y and connected_z):
print('not connected solids in all directions')
#-------------------------------------------------------------------------------
# Minkowski functionals
#-------------------------------------------------------------------------------
print("\nCompute the Minkowski functionals")
M0, M1, M2, M3 = compute_minkowski(M_bin)
print(f'M0 (porosity) = {M0:.3f}, M1 (specific surface area) = {M1:.3e}, M2 (mean grain size) = {M2:.3e}, M3 (Euler characteristic) = {M3:.3e} \n')
#-------------------------------------------------------------------------------
# fmm
#-------------------------------------------------------------------------------
# compute the geometrical tortuosities
print("Computing tortuosity on x for solid")
tau_x = compute_tortuosity_fast_marching(M_bin_extended_connected_x, extraction=[dim_sample, 2*dim_sample, dim_sample, 2*dim_sample], dx=1.0, neighborhood=6)
print("Computing tortuosity on y for solid")
tau_y = compute_tortuosity_fast_marching(np.transpose(M_bin_extended_connected_y, (1, 2, 0)), extraction=[dim_sample, 2*dim_sample, dim_sample, 2*dim_sample], dx=1.0, neighborhood=6)
print("Computing tortuosity on z for solid")
tau_z = compute_tortuosity_fast_marching(np.transpose(M_bin_extended_connected_z, (2, 0, 1)), extraction=[dim_sample, 2*dim_sample, dim_sample, 2*dim_sample], dx=1.0, neighborhood=6)
print(f"tau_x={tau_x:.3f}, tau_y={tau_y:.3f}, tau_z={tau_z:.3f}\n")
print("Computing tortuosity on x for pore")
tau_x = compute_tortuosity_fast_marching(M_inv_bin_extended_connected_x, extraction=[dim_sample, 2*dim_sample, dim_sample, 2*dim_sample], dx=1.0, neighborhood=6)
print("Computing tortuosity on y for pore")
tau_y = compute_tortuosity_fast_marching(np.transpose(M_inv_bin_extended_connected_y, (1, 2, 0)), extraction=[dim_sample, 2*dim_sample, dim_sample, 2*dim_sample], dx=1.0, neighborhood=6)
print("Computing tortuosity on z for pore")
tau_z = compute_tortuosity_fast_marching(np.transpose(M_inv_bin_extended_connected_z, (2, 0, 1)), extraction=[dim_sample, 2*dim_sample, dim_sample, 2*dim_sample], dx=1.0, neighborhood=6)
print(f"tau_x={tau_x:.3f}, tau_y={tau_y:.3f}, tau_z={tau_z:.3f}")