diff --git a/experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia.py b/experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia.py new file mode 100644 index 0000000..a205650 --- /dev/null +++ b/experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia.py @@ -0,0 +1,212 @@ +import sys +import os + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../..'))) + +from mpi4py import MPI +import numpy as np +import time +import matplotlib.pyplot as plt +from muGrid import Solvers + +from muFFTTO import domain +from muFFTTO import microstructure_library +from muFFTTO.visualization_utils import plot_field_on_grid + +from muFFTTO.grid_adaptation_methods import adapt_grid_to_circle, adapt_grid_to_circle_EXAMPLE_ + +# Copy of an example of how to usu muFFTTO to solve the homogenization problem for 2D heat conductivity problem +# using deformed grid with Jia's function for deformation + +problem_type = 'conductivity' +discretization_type = 'finite_element' +element_type = 'linear_triangles' +geometry_ID = 'square_inclusion' + +domain_size = (1, 1) +number_of_pixels = (32, 32) + +my_cell = domain.PeriodicUnitCell(domain_size=domain_size, + problem_type=problem_type) + +discretization = domain.Discretization(cell=my_cell, + nb_of_pixels_global=number_of_pixels, + discretization_type=discretization_type, + element_type=element_type) +start_time = time.time() + +# create material data field +mat_contrast = 1 +mat_contrast_2 = 1e2 +conductivity_C_1 = np.array([[1., 0], [0, 1.0]]) + +material_data_field_C_0 = discretization.get_material_data_size_field_mugrid(name='conductivity_tensor') + +# populate the field with C_1 material +material_data_field_C_0.s[...] = conductivity_C_1[:, :, np.newaxis, np.newaxis, np.newaxis] + +# material distribution +# TODO[Jia]: Here you have to place your function +result = adapt_grid_to_circle_EXAMPLE_( + nb_grid_points=number_of_pixels, + domain_size=domain_size, center=(0.5, 0.5), radius=0.2, + reference_grid_points_coords=discretization.get_nodal_points_coordinates().s[:, 0, ...], + iters=80, omega=0.8 +) + +coords_of_displaced_nodes = result["coords_of_displaced_nodes"] +phase_indicator_array = result["inside"] + +phase_field = discretization.get_scalar_field(name='phase_field') +phase_field.s[0, 0] = phase_indicator_array +matrix_mask = phase_indicator_array > 0 +inc_mask = phase_indicator_array == 0 + +# apply material distribution +material_data_field_C_0.s[..., matrix_mask] = mat_contrast_2 * material_data_field_C_0.s[..., matrix_mask] +material_data_field_C_0.s[..., inc_mask] = mat_contrast * material_data_field_C_0.s[..., inc_mask] +# --------------------------------------------------------------------------------------------------------------------- # +# Reference coordinates +ref_grid_coords_ixyz = discretization.get_nodal_points_coordinates().s[:, 0, ...] + +# Deformed coordinates +def_grid_coords_inxyz = discretization.get_displacement_sized_field(name='deformed_nodal_points_coordinates_inxyz') + +# grid_nodes_displacement +grid_nodes_displacement_inxyz = discretization.get_displacement_sized_field(name='grid_nodes_displacement_inxyz') + +grid_nodes_displacement_inxyz.s.fill(0) +grid_nodes_displacement_inxyz.s[:, 0, ...] = coords_of_displaced_nodes - ref_grid_coords_ixyz + +# fill in the deformation with analytical +def_grid_coords_inxyz.s[:, 0, ...] = ref_grid_coords_ixyz[...] + grid_nodes_displacement_inxyz.s[:, 0, ...] + +# Deformed coords with periodic extension for plotting +x_plot = discretization.get_nodal_points_coordinates_with_periodic_nodes() +# add deformation +x_plot[..., :-1, :-1] += grid_nodes_displacement_inxyz.s[...] +x_plot = np.squeeze(x_plot, axis=1) # removes axis for more nodal points +# Visualize grid and material +plot_field_on_grid(coordinates_for_plot=x_plot, field_to_plot=phase_field.s[0, 0], name='Material') + +# Deformation gradient F = I + grad(u) +F_ijqxy = discretization.get_displacement_gradient_sized_field(name='Grid_Deformation_gradient_F_ijqxy') +discretization.fft.communicate_ghosts(grid_nodes_displacement_inxyz) +discretization.apply_gradient_operator_mugrid(grid_nodes_displacement_inxyz, F_ijqxy) +F_ijqxy.s[...] += np.eye(2)[:, :, None, None, None] +# determinant and inverse of the deformation gradient +det_F = np.linalg.det(F_ijqxy.s.transpose(2, 3, 4, 0, 1)) +inv_F = np.linalg.pinv(F_ijqxy.s.transpose(2, 3, 4, 0, 1)).transpose(3, 4, 0, 1, 2) + +# plot def_F in grid +plot_field_on_grid(coordinates_for_plot=x_plot, field_to_plot=det_F[0], name='det(F)') + + +def K_fun(x, Ax): + """ + Matrix-free application of the Hessian matrix. For deformed grids + """ + + discretization.apply_system_matrix_mugrid_deformed_grid(material_data_field=material_data_field_C_0, + input_field_inxyz=x, + output_field_inxyz=Ax, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + discretization.fft.communicate_ghosts(Ax) + + +preconditioner = discretization.get_preconditioner_Green_mugrid(reference_material_data_ijkl=conductivity_C_1) + + +def M_fun(x, Px): + """ + Function to compute the product of the Preconditioner matrix with a vector. + The Preconditioner is represented by the convolution operator. + """ + discretization.fft.communicate_ghosts(x) + discretization.apply_preconditioner_mugrid(preconditioner_Fourier_fnfnqks=preconditioner, + input_nodal_field_fnxyz=x, + output_nodal_field_fnxyz=Px) + + +solution_field = discretization.get_unknown_size_field(name='solution') +macro_gradient_field = discretization.get_gradient_size_field(name='macro_gradient_field') +rhs_field = discretization.get_unknown_size_field(name='rhs_field') + +dim = discretization.domain_dimension +homogenized_A_ij = np.zeros(np.array(2 * [dim, ])) + +for i in range(dim): + # set macroscopic gradient + macro_gradient = np.zeros([dim]) + macro_gradient[i] = 1 + + macro_gradient_field.sg.fill(0) + discretization.get_macro_gradient_field_mugrid(macro_gradient_ij=macro_gradient, + macro_gradient_field_ijqxyz=macro_gradient_field) + + # Macro gradient in reference domain + macro_gradient_field.s[...] = np.einsum('ij...,jk...->ik...', macro_gradient_field.s[...], inv_F) + discretization.fft.communicate_ghosts(field=macro_gradient_field) + + # Solve equilibrium + rhs_field.sg.fill(0) + discretization.get_rhs_mugrid_deformed_grid(material_data_field_ijklqxyz=material_data_field_C_0, + macro_gradient_field_ijqxyz=macro_gradient_field, + rhs_inxyz=rhs_field, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + + + def callback(iteration, fields): + """ + Callback function to print the current solution, residual, and search direction. + """ + norm_of_rr = fields['rr'] + if discretization.communicator.rank == 0: + print(f"{iteration:5} norm of residual = {norm_of_rr:.5}") + + + Solvers.conjugate_gradients( + comm=discretization.communicator, + fc=discretization.field_collection, + hessp=K_fun, # linear operator + b=rhs_field, # right-hand side + x=solution_field, + prec=M_fun, + rtol=1e-6, + maxiter=2000, + callback=callback) + + if discretization.communicator.size == 1: + # Plot the first component of the solution field + plot_field_on_grid(coordinates_for_plot=x_plot, + field_to_plot=solution_field.s[0, 0], + name=f'Solution field - macro gradient {macro_gradient} ') + discretization.fft.communicate_ghosts(field=solution_field) + + sum_sol = discretization.mpi_reduction.sum(solution_field.s, + axis=tuple(range(-3, 0))) + print('rank' f'{MPI.COMM_WORLD.rank:6} sum_sol =' f'{sum_sol}') + + homogenized_A_ij[i, :] = discretization.get_homogenized_stress_mugrid_deformed_grid( + material_data_field_ijklqxyz=material_data_field_C_0, + temperature_field_inxyz=solution_field, + macro_gradient_field_ijqxyz=macro_gradient_field, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + + # ---------------------------------------------------------------------- + print( + "homogenized conductivity tangent =\n" + + np.array2string(homogenized_A_ij, formatter={'float_kind': lambda x: f"{x:0.8f}"}) + ) + +end_time = time.time() +elapsed_time = end_time - start_time +if discretization.communicator.rank == 0: + print("Elapsed time: ", elapsed_time, 'seconds') + print("Elapsed time: ", elapsed_time / 60, 'minutes') + J_eff = mat_contrast_2 * np.sqrt((mat_contrast_2 + 3 * mat_contrast) / (3 * mat_contrast_2 + mat_contrast)) + print(f'Analytical solution conductivity - A^eff_11 : {J_eff:0.8f}') + print(f'Numerical solution conductivity - A^eff_11 : {homogenized_A_ij[0, 0]:0.8f}') diff --git a/experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia_multiple_phases.py b/experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia_multiple_phases.py new file mode 100644 index 0000000..f2050a1 --- /dev/null +++ b/experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia_multiple_phases.py @@ -0,0 +1,217 @@ +import sys +import os + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../..'))) + +from mpi4py import MPI +import numpy as np +import time +import matplotlib.pyplot as plt +from muGrid import Solvers + +from muFFTTO import domain +from muFFTTO import microstructure_library +from muFFTTO.visualization_utils import plot_field_on_grid + +from muFFTTO.grid_adaptation_methods import adapt_grid_to_circle, adapt_grid_to_circle_EXAMPLE_ + +# Copy of an example of how to usu muFFTTO to solve the homogenization problem for 2D heat conductivity problem +# using deformed grid with Jia's function for deformation + +problem_type = 'conductivity' +discretization_type = 'finite_element' +element_type = 'linear_triangles' +geometry_ID = 'square_inclusion' + +domain_size = (1, 1) +number_of_pixels = (32, 32) + +my_cell = domain.PeriodicUnitCell(domain_size=domain_size, + problem_type=problem_type) + +discretization = domain.Discretization(cell=my_cell, + nb_of_pixels_global=number_of_pixels, + discretization_type=discretization_type, + element_type=element_type) +start_time = time.time() + +# create material data field +conductivity_C_0 = np.array([[1., 0], + [0, 1.0]]) +conductivity_C_1 = np.array([[2., 0], + [0, 2.0]]) +conductivity_C_2 = np.array([[3., 0], [0, 3.0]]) + +material_data_field_C_global = discretization.get_material_data_size_field_mugrid(name='conductivity_tensor') +material_data_field_C_global.s.fill(0) + +# populate the field with zeros material + +# material distribution +# TODO[Jia]: Here you have to place your function +result = adapt_grid_to_circle_EXAMPLE_( + nb_grid_points=number_of_pixels, + domain_size=domain_size, center=(0.5, 0.5), radius=0.2, + reference_grid_points_coords=discretization.get_nodal_points_coordinates().s[:, 0, ...], + iters=80, omega=0.8 +) + +coords_of_displaced_nodes = result["coords_of_displaced_nodes"] +phase_indicator_array = result["inside"] +# artificial adding of third phase +phase_indicator_array[5:10,7:16]=2 + +phase_field = discretization.get_scalar_field(name='phase_field') +phase_field.s[0, 0] = phase_indicator_array +# here you create masks for individual phases +mat_0_mask= phase_indicator_array == 0 +mat_1_mask= phase_indicator_array == 1 +mat_2_mask= phase_indicator_array == 2 + +# apply material distribution - assign specific local (d,d) conductivity matrix to the global data field +material_data_field_C_global.s[..., mat_0_mask] = conductivity_C_0[...,np.newaxis,np.newaxis] +material_data_field_C_global.s[..., mat_1_mask] = conductivity_C_1[...,np.newaxis,np.newaxis] +material_data_field_C_global.s[..., mat_2_mask] = conductivity_C_2[...,np.newaxis,np.newaxis] +# --------------------------------------------------------------------------------------------------------------------- # +# Reference coordinates +ref_grid_coords_ixyz = discretization.get_nodal_points_coordinates().s[:, 0, ...] + +# Deformed coordinates +def_grid_coords_inxyz = discretization.get_displacement_sized_field(name='deformed_nodal_points_coordinates_inxyz') + +# grid_nodes_displacement +grid_nodes_displacement_inxyz = discretization.get_displacement_sized_field(name='grid_nodes_displacement_inxyz') + +grid_nodes_displacement_inxyz.s.fill(0) +grid_nodes_displacement_inxyz.s[:, 0, ...] = coords_of_displaced_nodes - ref_grid_coords_ixyz + +# fill in the deformation with analytical +def_grid_coords_inxyz.s[:, 0, ...] = ref_grid_coords_ixyz[...] + grid_nodes_displacement_inxyz.s[:, 0, ...] + +# Deformed coords with periodic extension for plotting +x_plot = discretization.get_nodal_points_coordinates_with_periodic_nodes() +# add deformation +x_plot[..., :-1, :-1] += grid_nodes_displacement_inxyz.s[...] +x_plot = np.squeeze(x_plot, axis=1) # removes axis for more nodal points +# Visualize grid and material +plot_field_on_grid(coordinates_for_plot=x_plot, field_to_plot=phase_field.s[0, 0], name='Material') + +# Deformation gradient F = I + grad(u) +F_ijqxy = discretization.get_displacement_gradient_sized_field(name='Grid_Deformation_gradient_F_ijqxy') +discretization.fft.communicate_ghosts(grid_nodes_displacement_inxyz) +discretization.apply_gradient_operator_mugrid(grid_nodes_displacement_inxyz, F_ijqxy) +F_ijqxy.s[...] += np.eye(2)[:, :, None, None, None] +# determinant and inverse of the deformation gradient +det_F = np.linalg.det(F_ijqxy.s.transpose(2, 3, 4, 0, 1)) +inv_F = np.linalg.pinv(F_ijqxy.s.transpose(2, 3, 4, 0, 1)).transpose(3, 4, 0, 1, 2) + +# plot def_F in grid +plot_field_on_grid(coordinates_for_plot=x_plot, field_to_plot=det_F[0], name='det(F)') + + +def K_fun(x, Ax): + """ + Matrix-free application of the Hessian matrix. For deformed grids + """ + + discretization.apply_system_matrix_mugrid_deformed_grid(material_data_field=material_data_field_C_global, + input_field_inxyz=x, + output_field_inxyz=Ax, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + discretization.fft.communicate_ghosts(Ax) + + +preconditioner = discretization.get_preconditioner_Green_mugrid(reference_material_data_ijkl=conductivity_C_1) + + +def M_fun(x, Px): + """ + Function to compute the product of the Preconditioner matrix with a vector. + The Preconditioner is represented by the convolution operator. + """ + discretization.fft.communicate_ghosts(x) + discretization.apply_preconditioner_mugrid(preconditioner_Fourier_fnfnqks=preconditioner, + input_nodal_field_fnxyz=x, + output_nodal_field_fnxyz=Px) + + +solution_field = discretization.get_unknown_size_field(name='solution') +macro_gradient_field = discretization.get_gradient_size_field(name='macro_gradient_field') +rhs_field = discretization.get_unknown_size_field(name='rhs_field') + +dim = discretization.domain_dimension +homogenized_A_ij = np.zeros(np.array(2 * [dim, ])) + +for i in range(dim): + # set macroscopic gradient + macro_gradient = np.zeros([dim]) + macro_gradient[i] = 1 + + macro_gradient_field.sg.fill(0) + discretization.get_macro_gradient_field_mugrid(macro_gradient_ij=macro_gradient, + macro_gradient_field_ijqxyz=macro_gradient_field) + + # Macro gradient in reference domain + macro_gradient_field.s[...] = np.einsum('ij...,jk...->ik...', macro_gradient_field.s[...], inv_F) + discretization.fft.communicate_ghosts(field=macro_gradient_field) + + # Solve equilibrium + rhs_field.sg.fill(0) + discretization.get_rhs_mugrid_deformed_grid(material_data_field_ijklqxyz=material_data_field_C_global, + macro_gradient_field_ijqxyz=macro_gradient_field, + rhs_inxyz=rhs_field, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + + + def callback(iteration, fields): + """ + Callback function to print the current solution, residual, and search direction. + """ + norm_of_rr = fields['rr'] + if discretization.communicator.rank == 0: + print(f"{iteration:5} norm of residual = {norm_of_rr:.5}") + + + Solvers.conjugate_gradients( + comm=discretization.communicator, + fc=discretization.field_collection, + hessp=K_fun, # linear operator + b=rhs_field, # right-hand side + x=solution_field, + prec=M_fun, + rtol=1e-6, + maxiter=2000, + callback=callback) + + if discretization.communicator.size == 1: + # Plot the first component of the solution field + plot_field_on_grid(coordinates_for_plot=x_plot, + field_to_plot=solution_field.s[0, 0], + name=f'Solution field - macro gradient {macro_gradient} ') + discretization.fft.communicate_ghosts(field=solution_field) + + sum_sol = discretization.mpi_reduction.sum(solution_field.s, + axis=tuple(range(-3, 0))) + print('rank' f'{MPI.COMM_WORLD.rank:6} sum_sol =' f'{sum_sol}') + + homogenized_A_ij[i, :] = discretization.get_homogenized_stress_mugrid_deformed_grid( + material_data_field_ijklqxyz=material_data_field_C_global, + temperature_field_inxyz=solution_field, + macro_gradient_field_ijqxyz=macro_gradient_field, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + + # ---------------------------------------------------------------------- + print( + "homogenized conductivity tangent =\n" + + np.array2string(homogenized_A_ij, formatter={'float_kind': lambda x: f"{x:0.8f}"}) + ) + +end_time = time.time() +elapsed_time = end_time - start_time +if discretization.communicator.rank == 0: + print("Elapsed time: ", elapsed_time, 'seconds') + print("Elapsed time: ", elapsed_time / 60, 'minutes') + print(f'Numerical solution conductivity - A^eff_11 : {homogenized_A_ij[0, 0]:0.8f}') diff --git a/experiments/grid_adaptation/exp_2D_homogenization_elasticity_transformed_grid_Jia.py b/experiments/grid_adaptation/exp_2D_homogenization_elasticity_transformed_grid_Jia.py new file mode 100644 index 0000000..83a7b10 --- /dev/null +++ b/experiments/grid_adaptation/exp_2D_homogenization_elasticity_transformed_grid_Jia.py @@ -0,0 +1,233 @@ +import sys +import os + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../..'))) + +from mpi4py import MPI +import numpy as np +import time +import matplotlib.pyplot as plt +from muGrid import Solvers + +from muFFTTO import domain +from muFFTTO import microstructure_library + +from muFFTTO.visualization_utils import plot_field_on_grid, get_deformed_grid_coords_two_dim +from muFFTTO.grid_adaptation_methods import adapt_grid_to_circle, adapt_grid_to_circle_EXAMPLE_ + +# Example of how to usu muFFTTO to solve the homogenization problem for 2D elasticity problem +# using deformed grid + +problem_type = 'elasticity' +discretization_type = 'finite_element' +element_type = 'linear_triangles' +formulation = 'small_strain' +geometry_ID = 'square_inclusion' + +domain_size = (1, 1) +number_of_pixels = (32, 32) + +my_cell = domain.PeriodicUnitCell(domain_size=domain_size, + problem_type=problem_type) + +discretization = domain.Discretization(cell=my_cell, + nb_of_pixels_global=number_of_pixels, + discretization_type=discretization_type, + element_type=element_type) +start_time = time.time() +# initialize material data +K_0, G_0 = domain.get_bulk_and_shear_modulus(E=1, poison=0.2) + +# create material data field +elastic_C_1 = domain.get_elastic_material_tensor(dim=discretization.domain_dimension, + K=K_0, + mu=G_0, + kind='linear') +if discretization.communicator.rank == 0: + print('elastic tangent = \n {}'.format(domain.compute_Voigt_notation_4order(elastic_C_1))) + +material_data_field_C_0 = discretization.get_material_data_size_field_mugrid(name='elastic_tensor') + +# populate the field with C_1 material +material_data_field_C_0.s[...] = elastic_C_1[:, :, :, :, np.newaxis, np.newaxis, np.newaxis] + +# material distribution +phase_field = discretization.get_scalar_field(name='phase_field') + +# TODO[Jia]: Here you have to place your function +result = adapt_grid_to_circle_EXAMPLE_( + nb_grid_points=number_of_pixels, + domain_size=domain_size, center=(0.5, 0.5), radius=0.2, + reference_grid_points_coords=discretization.get_nodal_points_coordinates().s[:, 0, ...], + iters=80, omega=0.8 +) + +coords_of_displaced_nodes = result["coords_of_displaced_nodes"] +phase_indicator_array = result["inside"] + + +phase_field.s[0, 0] = phase_indicator_array + +mat_contrast = 1e2 +mat_contrast_2 = 1 +matrix_mask = phase_field.s[0, 0] > 0 +inc_mask = phase_field.s[0, 0] == 0 + +# apply material distribution +material_data_field_C_0.s[..., matrix_mask] = mat_contrast_2 * material_data_field_C_0.s[..., matrix_mask] +material_data_field_C_0.s[..., inc_mask] = mat_contrast * material_data_field_C_0.s[..., inc_mask] + +# --------------------------------------------------------------------------------------------------------------------- # +# Reference coordinates +ref_grid_coords_ixyz = discretization.fft.coords + +# Deformed coordinates +def_grid_coords_inxyz = discretization.get_displacement_sized_field(name='deformed_nodal_points_coordinates_inxyz') + +# grid_nodes_displacement +grid_nodes_displacement_inxyz = discretization.get_displacement_sized_field(name='grid_nodes_displacement_inxyz') + +grid_nodes_displacement_inxyz.s.fill(0) +grid_nodes_displacement_inxyz.s[:, 0, ...] = coords_of_displaced_nodes - ref_grid_coords_ixyz + +# fill in the deformation with analytical +def_grid_coords_inxyz.s[:, 0, ...] = ref_grid_coords_ixyz[...] + grid_nodes_displacement_inxyz.s[:, 0, ...] + +# Deformed coords with periodic extension for plotting +x_plot = discretization.get_nodal_points_coordinates_with_periodic_nodes() +# add deformation +x_plot[..., :-1, :-1] += grid_nodes_displacement_inxyz.s[...] +x_plot = np.squeeze(x_plot, axis=1) # removes axis for more nodal points +# Visualize grid and material +plot_field_on_grid(coordinates_for_plot=x_plot, field_to_plot=phase_field.s[0, 0], name='Material') + +# Deformation gradient F = I + grad(u) +F_ijqxy = discretization.get_displacement_gradient_sized_field(name='Grid_Deformation_gradient_F_ijqxy') +discretization.fft.communicate_ghosts(grid_nodes_displacement_inxyz) +discretization.apply_gradient_operator_mugrid(grid_nodes_displacement_inxyz, F_ijqxy) +F_ijqxy.s[...] += np.eye(2)[:, :, None, None, None] +# determinant and inverse of the deformation gradient +det_F = np.linalg.det(F_ijqxy.s.transpose(2, 3, 4, 0, 1)) +inv_F = np.linalg.pinv(F_ijqxy.s.transpose(2, 3, 4, 0, 1)).transpose(3, 4, 0, 1, 2) + +# plot def_F in grid +plot_field_on_grid(coordinates_for_plot=x_plot, field_to_plot=det_F[0], name='det(F)') + + +# -------------------------------------------------- + +def K_fun(x, Ax): + """ + Matrix-free application of the Hessian matrix. For deformed grids + """ + + discretization.apply_system_matrix_mugrid_deformed_grid(material_data_field=material_data_field_C_0, + input_field_inxyz=x, + output_field_inxyz=Ax, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F, + formulation='small_strain') + discretization.fft.communicate_ghosts(Ax) + + +preconditioner = discretization.get_preconditioner_Green_mugrid(reference_material_data_ijkl=elastic_C_1) + + +def M_fun(x, Px): + """ + Function to compute the product of the Preconditioner matrix with a vector. + The Preconditioner is represented by the convolution operator. + """ + discretization.fft.communicate_ghosts(x) + discretization.apply_preconditioner_mugrid(preconditioner_Fourier_fnfnqks=preconditioner, + input_nodal_field_fnxyz=x, + output_nodal_field_fnxyz=Px) + + +displacement_fluctuation = discretization.get_unknown_size_field(name='displacement_fluctuation') +macro_gradient_field = discretization.get_gradient_size_field(name='macro_gradient_field') +rhs_field = discretization.get_unknown_size_field(name='rhs_field') + +dim = discretization.domain_dimension +homogenized_C_ijkl = np.zeros(np.array(4 * [dim, ])) +# compute whole homogenized elastic tangent +for i in range(dim): + for j in range(dim): + # set macroscopic gradient + macro_gradient_ij = np.zeros([dim, dim]) + macro_gradient_ij[i, j] = 0.3 + + macro_gradient_field.sg.fill(0) + discretization.get_macro_gradient_field_mugrid(macro_gradient_ij=macro_gradient_ij, + macro_gradient_field_ijqxyz=macro_gradient_field) + + # Macro gradient in reference domain + macro_gradient_field.s[...] = np.einsum('ij...,jk...->ik...', macro_gradient_field.s[...], inv_F) + discretization.fft.communicate_ghosts(field=macro_gradient_field) + + # Solve mechanical equilibrium constrain + rhs_field.sg.fill(0) + discretization.get_rhs_mugrid_deformed_grid(material_data_field_ijklqxyz=material_data_field_C_0, + macro_gradient_field_ijqxyz=macro_gradient_field, + rhs_inxyz=rhs_field, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + + + def callback(iteration, fields): + """ + Callback function to print the current solution, residual, and search direction. + """ + norm_of_rr = fields['rr'] + if discretization.communicator.rank == 0: + print(f"{iteration:5} norm of residual = {norm_of_rr:.5}") + + + Solvers.conjugate_gradients( + comm=discretization.communicator, + fc=discretization.field_collection, + hessp=K_fun, # linear operator + b=rhs_field, # right-hand side + x=displacement_fluctuation, + prec=M_fun, + tol=1e-6, + maxiter=2000, + callback=callback) + + if discretization.communicator.size == 1: + # plot deformed domain + x_plot_ixyz = get_deformed_grid_coords_two_dim(discretization=discretization, + grid_nodes_displacement_inxyz=grid_nodes_displacement_inxyz, + macro_gradient_ij=macro_gradient_ij, + displacement_fluctuation=displacement_fluctuation) + + x_plot_ixyz[1, -1, -1] += displacement_fluctuation.s[1, 0, 0, 0] + + # ref_grid_coords_ixyz + plot_field_on_grid(coordinates_for_plot=x_plot_ixyz, + field_to_plot=displacement_fluctuation.s[0, 0], + name=f'Solution field - macro gradient {macro_gradient_ij} ') + + discretization.fft.communicate_ghosts(field=displacement_fluctuation) + + sum_sol = discretization.mpi_reduction.sum(displacement_fluctuation.s, + axis=tuple(range(-3, 0))) + print('rank' f'{MPI.COMM_WORLD.rank:6} sum_sol =' f'{sum_sol}') + + homogenized_C_ijkl[i, j] = discretization.get_homogenized_stress_mugrid_deformed_grid( + material_data_field_ijklqxyz=material_data_field_C_0, + temperature_field_inxyz=displacement_fluctuation, + macro_gradient_field_ijqxyz=macro_gradient_field, + det_of_deformation_gradient=det_F, + inv_of_deformation_gradient=inv_F) + + # ---------------------------------------------------------------------- + print( + "Homogenized elastic tangent =\n" + + np.array2string(domain.compute_Voigt_notation_4order(homogenized_C_ijkl), + formatter={'float_kind': lambda x: f"{x:0.8f}"}) + ) + +end_time = time.time() +elapsed_time = end_time - start_time +print("Elapsed time: ", elapsed_time) diff --git a/experiments/grid_adaptation/exp_analytical_grid_adaptation.py b/experiments/grid_adaptation/exp_analytical_grid_adaptation.py new file mode 100644 index 0000000..00bb3e0 --- /dev/null +++ b/experiments/grid_adaptation/exp_analytical_grid_adaptation.py @@ -0,0 +1,166 @@ +import argparse +import time +import numpy as np +import matplotlib.pyplot as plt +from muFFTTO.analytical_grid_adaptation import adapt_grid_to_circle + +# Fixed parameters +DOMAIN_HALF_SIZE = 40.0 +CIRCLE_CENTER = (0.0, 0.0) +CIRCLE_RADIUS = 20.0 +NB_ITERS = 400 +OMEGA = 0.8 +K0 = 0.25 +KMIN = 0.005 + +#coor = make_grid_nodes(N=8, L=2.0) +#print(coor) + +def parse_args() -> argparse.Namespace: + """ + Parse command-line arguments for the grid-adaptation experiment. + + The function reads optional command-line arguments specifying the mesh + size and the stiffness-decay exponent. If no values are provided in the + terminal, default values are used. + + Parameters + ---------- + None + + Returns + ------- + argparse.Namespace + Namespace containing the parsed command-line arguments. + - args.N : int + Mesh size used for the grid adaptation. + - args.b : float + Stiffness-decay exponent used for the grid adaptation. + """ + parser = argparse.ArgumentParser( + description="Run analytical grid adaptation for one specified N and b." + ) + parser.add_argument( + "--N", + type=int, + default=16, + help="Grid size N (default: 16)", + ) + parser.add_argument( + "--b", + type=float, + default=0.5, + help="Stiffness exponent b (default: 0.5)", + ) + return parser.parse_args() + +def plot_initial_and_adapted_mesh( + P0: np.ndarray, + P: np.ndarray, + N: int, + b: float, + elapsed_time: float, +) -> None: + """ + Plot the initial and adapted meshes side by side. + + The function visualizes the grid lines of the initial regular mesh and the + adapted mesh obtained after grid deformation toward a circular geometry. + The selected mesh size, stiffness exponent, and total runtime are shown in + the figure title. + + Parameters + ---------- + P0 : numpy.ndarray + Initial nodal coordinates with shape [nx, ny, xy]. + P : numpy.ndarray + Adapted nodal coordinates with shape [nx, ny, xy]. + N : int + Mesh size used for the grid adaptation. + b : float + Stiffness-decay exponent used for the grid adaptation. + elapsed_time : float + Total runtime of the grid-adaptation procedure in seconds. + + Returns + ------- + None + """ + N_plot = P0.shape[0] - 1 + + fig, axes = plt.subplots(1, 2, figsize=(11, 5)) + ax_init, ax_def = axes + + # Initial grid + for i in range(N_plot + 1): + ax_init.plot(P0[i, :, 0], P0[i, :, 1], color="0.2", linewidth=0.7) + for j in range(N_plot + 1): + ax_init.plot(P0[:, j, 0], P0[:, j, 1], color="0.2", linewidth=0.7) + + ax_init.set_aspect("equal", adjustable="box") + ax_init.set_title("Initial grid") + ax_init.set_xlabel("x") + ax_init.set_ylabel("y") + + # Adapted grid + for i in range(N_plot + 1): + ax_def.plot(P[i, :, 0], P[i, :, 1], color="0.2", linewidth=0.7) + for j in range(N_plot + 1): + ax_def.plot(P[:, j, 0], P[:, j, 1], color="0.2", linewidth=0.7) + + ax_def.set_aspect("equal", adjustable="box") + ax_def.set_title("Adapted grid") + ax_def.set_xlabel("x") + ax_def.set_ylabel("y") + + fig.suptitle( + f"Grid adaptation to circle | N = {N}, b = {b:.2f}, time = {elapsed_time:.4f} s" + ) + fig.tight_layout() + plt.show() + + +def main() -> None: + """ + Run one grid-adaptation experiment and display the resulting meshes. + + The function parses the command-line arguments, executes the analytical + grid-adaptation procedure for one pair of parameters ``N`` and ``b``, + measures the runtime, and visualizes the initial and adapted meshes. + + Parameters + ---------- + None + + Returns + ------- + None + """ + args = parse_args() + + start_time = time.perf_counter() + result = adapt_grid_to_circle( + N=args.N, + L=DOMAIN_HALF_SIZE, + center=CIRCLE_CENTER, + R=CIRCLE_RADIUS, + iters=NB_ITERS, + omega=OMEGA, + b=args.b, + k0=K0, + kmin=KMIN, + ) + end_time = time.perf_counter() + + elapsed_time = end_time - start_time + + plot_initial_and_adapted_mesh( + result["P0"], + result["P"], + N=args.N, + b=args.b, + elapsed_time=elapsed_time, + ) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/experiments/grid_adaptation/exp_circle_muGrid_real_field.py b/experiments/grid_adaptation/exp_circle_muGrid_real_field.py new file mode 100644 index 0000000..6c5752e --- /dev/null +++ b/experiments/grid_adaptation/exp_circle_muGrid_real_field.py @@ -0,0 +1,340 @@ +import argparse +import time +import numpy as np +import matplotlib.pyplot as plt + +from muFFTTO.grid_adaptation_methods import ( + adapt_grid_to_circle, + pack_adapted_grid_to_fields, + print_field_summary, +) + + +DOMAIN_SIZE = 40.0 +CIRCLE_CENTER = (20.0, 20.0) +CIRCLE_RADIUS = 10.0 +NB_ITERS = 400 +OMEGA = 0.8 +K0 = 0.25 +KMIN = 0.005 + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Run periodic-style analytical grid adaptation for one specified N and b." + ) + parser.add_argument( + "--N", + type=int, + default=32, + help="Number of cells / stored points per direction (default: 32)", + ) + parser.add_argument( + "--b", + type=float, + default=1.0, + help="Stiffness exponent b (default: 1.0)", + ) + parser.add_argument( + "--ghosts", + type=int, + default=1, + help="Number of muGrid ghost cells per side (default: 1)", + ) + parser.add_argument( + "--no-store-positions", + action="store_true", + help="Do not store P0 and P as muGrid fields", + ) + return parser.parse_args() + + +def close_periodic_grid_for_plot(P: np.ndarray, L: float) -> np.ndarray: + """ + Create a temporary closed grid for plotting. + + Input shape is [xy, N, N]. Output shape is [xy, N+1, N+1]. + The last column copies the first column with x shifted by +L. + The last row copies the first row with y shifted by +L. + """ + N = P.shape[1] + Pc = np.zeros((2, N + 1, N + 1), dtype=P.dtype) + + Pc[:, :N, :N] = P + Pc[:, N, :N] = P[:, 0, :] + Pc[0, N, :N] += L + + Pc[:, :N, N] = P[:, :, 0] + Pc[1, :N, N] += L + + Pc[:, N, N] = P[:, 0, 0] + Pc[0, N, N] += L + Pc[1, N, N] += L + + return Pc + + +def plot_initial_and_adapted_mesh( + P0: np.ndarray, + P: np.ndarray, + N: int, + b: float, + elapsed_time: float, + L: float, +) -> None: + P0c = close_periodic_grid_for_plot(P0, L=L) + Pc = close_periodic_grid_for_plot(P, L=L) + npx, npy = P0c.shape[1], P0c.shape[2] + + fig, axes = plt.subplots(1, 2, figsize=(11, 5)) + ax_init, ax_def = axes + + for i in range(npx): + ax_init.plot(P0c[0, i, :], P0c[1, i, :], color="0.2", linewidth=0.7) + for j in range(npy): + ax_init.plot(P0c[0, :, j], P0c[1, :, j], color="0.2", linewidth=0.7) + + ax_init.set_aspect("equal", adjustable="box") + ax_init.set_title("Initial periodic grid") + ax_init.set_xlabel("x") + ax_init.set_ylabel("y") + + for i in range(npx): + ax_def.plot(Pc[0, i, :], Pc[1, i, :], color="0.2", linewidth=0.7) + for j in range(npy): + ax_def.plot(Pc[0, :, j], Pc[1, :, j], color="0.2", linewidth=0.7) + + ax_def.set_aspect("equal", adjustable="box") + ax_def.set_title("Adapted periodic grid") + ax_def.set_xlabel("x") + ax_def.set_ylabel("y") + + fig.suptitle( + f"Periodic grid adaptation to circle | N = {N}, b = {b:.2f}, time = {elapsed_time:.4f} s" + ) + fig.tight_layout() + plt.show() + + +def print_numpy_summary(field_bundle: dict) -> None: + disp = field_bundle["numpy"]["displacement"] + inside = field_bundle["numpy"]["inside"] + interface = field_bundle["numpy"]["interface"] + dist = field_bundle["numpy"]["dist"] + k_node = field_bundle["numpy"]["k_node"] + + print("\nNumpy summary") + print("-" * 72) + print("displacement shape:", disp.shape) + print("ux min/max:", disp[0].min(), disp[0].max()) + print("uy min/max:", disp[1].min(), disp[1].max()) + print("inside sum:", inside.sum()) + print("interface sum:", interface.sum()) + print("dist min/max:", dist.min(), dist.max()) + print("k_node min/max:", k_node.min(), k_node.max()) + + +def main() -> None: + args = parse_args() + + start_time = time.perf_counter() + result = adapt_grid_to_circle( + N=args.N, + L=DOMAIN_SIZE, + center=CIRCLE_CENTER, + R=CIRCLE_RADIUS, + iters=NB_ITERS, + omega=OMEGA, + b=args.b, + k0=K0, + kmin=KMIN, + ) + end_time = time.perf_counter() + elapsed_time = end_time - start_time + + plot_initial_and_adapted_mesh( + result["P0"], + result["P"], + N=args.N, + b=args.b, + elapsed_time=elapsed_time, + L=DOMAIN_SIZE, + ) + + field_bundle = pack_adapted_grid_to_fields( + result, + ghosts=args.ghosts, + store_positions=not args.no_store_positions, + verbose=True, + ) + + print_field_summary(field_bundle) + print_numpy_summary(field_bundle) + #print(field_bundle) + + +if __name__ == "__main__": + main() + +''' +//successful version + +import argparse +import time +import numpy as np +import matplotlib.pyplot as plt + + +from muFFTTO.feature import adapt_grid_to_circle + + + +DOMAIN_SIZE = 40.0 +CIRCLE_CENTER = (20.0, 20.0) +CIRCLE_RADIUS = 10.0 +NB_ITERS = 400 +OMEGA = 0.8 +K0 = 0.25 +KMIN = 0.005 + + + +def parse_args() -> argparse.Namespace: +    parser = argparse.ArgumentParser( +        description="Run periodic-style analytical grid adaptation for one specified N and b." +    ) +    parser.add_argument( +        "--N", +        type=int, +        default=32, +        help="Number of cells / stored points per direction (default: 32)", +    ) +    parser.add_argument( +        "--b", +        type=float, +        default=1.0, +        help="Stiffness exponent b (default: 1.0)", +    ) +    return parser.parse_args() + + + + +def close_periodic_grid_for_plot(P: np.ndarray, L: float) -> np.ndarray: +    """ +    Create a temporary closed grid for plotting. + + +    Input shape is [xy, N, N]. Output shape is [xy, N+1, N+1]. +    The last column copies the first column with x shifted by +L. +    The last row copies the first row with y shifted by +L. +    """ +    N = P.shape[1] +    Pc = np.zeros((2, N + 1, N + 1), dtype=P.dtype) + + +    Pc[:, :N, :N] = P +    Pc[:, N, :N] = P[:, 0, :] +    Pc[0, N, :N] += L + + +    Pc[:, :N, N] = P[:, :, 0] +    Pc[1, :N, N] += L + + +    Pc[:, N, N] = P[:, 0, 0] +    Pc[0, N, N] += L +    Pc[1, N, N] += L + + +    return Pc + + + + +def plot_initial_and_adapted_mesh( +    P0: np.ndarray, +    P: np.ndarray, +    N: int, +    b: float, +    elapsed_time: float, +    L: float, +) -> None: +    P0c = close_periodic_grid_for_plot(P0, L=L) +    Pc = close_periodic_grid_for_plot(P, L=L) +    npx, npy = P0c.shape[1], P0c.shape[2] + + +    fig, axes = plt.subplots(1, 2, figsize=(11, 5)) +    ax_init, ax_def = axes + + +    for i in range(npx): +        ax_init.plot(P0c[0, i, :], P0c[1, i, :], color="0.2", linewidth=0.7) +    for j in range(npy): +        ax_init.plot(P0c[0, :, j], P0c[1, :, j], color="0.2", linewidth=0.7) + + +    ax_init.set_aspect("equal", adjustable="box") +    ax_init.set_title("Initial periodic grid") +    ax_init.set_xlabel("x") +    ax_init.set_ylabel("y") + + +    for i in range(npx): +        ax_def.plot(Pc[0, i, :], Pc[1, i, :], color="0.2", linewidth=0.7) +    for j in range(npy): +        ax_def.plot(Pc[0, :, j], Pc[1, :, j], color="0.2", linewidth=0.7) + + +    ax_def.set_aspect("equal", adjustable="box") +    ax_def.set_title("Adapted periodic grid") +    ax_def.set_xlabel("x") +    ax_def.set_ylabel("y") + + +    fig.suptitle( +        f"Periodic grid adaptation to circle | N = {N}, b = {b:.2f}, time = {elapsed_time:.4f} s" +    ) +    fig.tight_layout() +    plt.show() + + + + +def main() -> None: +    args = parse_args() + + +    start_time = time.perf_counter() +    result = adapt_grid_to_circle( +        N=args.N, +        L=DOMAIN_SIZE, +        center=CIRCLE_CENTER, +        R=CIRCLE_RADIUS, +        iters=NB_ITERS, +        omega=OMEGA, +        b=args.b, +        k0=K0, +        kmin=KMIN, +    ) +    end_time = time.perf_counter() + + +    elapsed_time = end_time - start_time + + +    plot_initial_and_adapted_mesh( +        result["P0"], +        result["P"], +        N=args.N, +        b=args.b, +        elapsed_time=elapsed_time, +        L=DOMAIN_SIZE, +    ) + + + +if __name__ == "__main__": +    main() +''' \ No newline at end of file diff --git a/muFFTTO/otsu.py b/muFFTTO/otsu.py new file mode 100644 index 0000000..d14be0c --- /dev/null +++ b/muFFTTO/otsu.py @@ -0,0 +1,421 @@ +from __future__ import annotations + +from typing import Any + +import cv2 +import numpy as np + + +def _normalize_to_u8(data: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Normalize a 2D image to [0, 1] and uint8 [0, 255].""" + data = np.asarray(data, dtype=np.float32) + + if data.ndim != 2: + raise ValueError(f"Expected a 2D image, got shape {data.shape}.") + + data_min = float(data.min()) + data_max = float(data.max()) + value_range = data_max - data_min + + if value_range <= np.finfo(np.float32).eps: + image_norm = np.zeros_like(data, dtype=np.float32) + else: + image_norm = (data - data_min) / value_range + + image_u8 = np.clip(image_norm * 255.0, 0, 255).astype(np.uint8) + return image_norm, image_u8 + + +def _validate_odd_positive(name: str, value: int) -> None: + """Validate a positive odd OpenCV kernel size.""" + if value <= 0 or value % 2 == 0: + raise ValueError(f"{name} must be a positive odd integer, got {value}.") + + +def otsu_edgeDetection_and_phaseIndicator( + data: np.ndarray, + blur_ksize: int = 7, + blur_sigma: float = 1.5, + morph_kernel_size: int = 3, + open_iterations: int = 1, + close_iterations: int = 1, + regions_label: bool = True, + connectivity: int = 8, + return_intermediate: bool = False, + +) -> dict[str, Any]: + """ + Segment dark regions in a two-dimensional grayscale image, extract their + external boundaries, and optionally assign a unique integer label to each + disconnected foreground region. + + The input image is normalized to uint8 intensity values, smoothed with a + Gaussian filter, and segmented using inverted Otsu thresholding. Morphological + opening and closing are then applied to reduce isolated noise, remove small + foreground artifacts, fill small holes, and close short gaps. + + Inverted Otsu thresholding treats dark pixels as foreground: + + - ``phase_mask_binary == 1`` identifies dark segmented regions. + - ``phase_mask_binary == 0`` identifies bright/background regions. + + When ``regions_label=True``, connected-component analysis assigns one unique + integer ID to each disconnected foreground region. Background pixels always + have label ``0``; foreground region labels are consecutive integers + ``1, 2, ..., N``. + + This function currently distinguishes regions by geometric connectivity only. + Therefore, two disconnected regions receive different labels even when their + average grayscale values are similar. A later processing step may group such + regions into shared material classes based on their mean grayscale values. + + Parameters + ---------- + data : numpy.ndarray + Two-dimensional grayscale input image with shape ``(nx, ny)``. Integer and + floating-point input dtypes are supported. + + blur_ksize : int, default=7 + Positive odd Gaussian-kernel width and height. Larger values produce + stronger smoothing before Otsu thresholding. + + blur_sigma : float, default=1.5 + Standard deviation of the Gaussian filter. Must be non-negative. + + morph_kernel_size : int, default=3 + Positive odd side length of the square structuring element used for + morphological opening and closing. + + open_iterations : int, default=1 + Number of morphological-opening iterations. Opening suppresses small + isolated foreground pixels and thin noise. + + close_iterations : int, default=1 + Number of morphological-closing iterations. Closing fills small holes and + connects short gaps inside foreground regions. + + regions_label : bool, default=True + If ``True``, perform connected-component labeling on the cleaned binary + foreground mask. + + - ``True``: ``phase_mask_label`` has dtype ``int32`` and contains labels + ``0, 1, ..., N``. + - ``False``: ``phase_mask_label`` is a uint8 copy of + ``phase_mask_binary`` and contains only ``0`` and ``1``. + + connectivity : int, default=8 + Connectivity used for connected-component labeling when + ``regions_label=True``. + + - ``4`` connects pixels sharing a horizontal or vertical edge. + - ``8`` additionally connects pixels touching diagonally at a corner. + + return_intermediate : bool, default=False + If ``True``, include normalized, blurred, thresholded, and morphology + intermediate arrays in the returned dictionary. + + Returns + ------- + dict[str, Any] + Dictionary containing the following arrays with shape ``(nx, ny)``: + + ``edge_mask`` : numpy.ndarray, uint8 + One-pixel-wide mask of external contours. A value of ``1`` marks an + outer boundary of a segmented foreground region. + + ``phase_mask_binary`` : numpy.ndarray, uint8 + Cleaned binary segmentation mask with values ``0`` for background and + ``1`` for dark foreground regions. + + ``phase_mask_label`` : numpy.ndarray + Region-label field. + + - If ``regions_label=True``, dtype is ``int32`` and values are + ``0, 1, ..., N``. + - If ``regions_label=False``, dtype is ``uint8`` and values are + ``0`` and ``1``. + + When ``regions_label=True``, the dictionary additionally contains: + + ``number_of_phase_regions`` : int + Number of disconnected foreground regions. Background label ``0`` is + excluded. + + ``phase_region_stats`` : numpy.ndarray + OpenCV connected-component statistics. Each row corresponds to one + region label, including background label ``0``. The columns contain + bounding-box x/y coordinates, width, height, and pixel area. + + ``phase_region_centroids`` : numpy.ndarray + Centroid coordinates ``(x, y)`` for every label, including background. + + When ``return_intermediate=True``, the dictionary additionally contains: + + ``image_norm`` : numpy.ndarray + Input image normalized to the interval ``[0, 1]``. + + ``image_u8`` : numpy.ndarray + Normalized uint8 image with intensity values in ``[0, 255]``. + + ``image_blur`` : numpy.ndarray + Gaussian-smoothed uint8 image. + + ``mask_raw`` : numpy.ndarray, uint8 + Binary result directly after inverted Otsu thresholding. + + ``mask_open`` : numpy.ndarray, uint8 + Binary result after morphological opening. + + ``mask_clean`` : numpy.ndarray, uint8 + Binary result after opening and closing; equal to + ``phase_mask_binary``. + + ``otsu_threshold`` : float + Automatically selected Otsu threshold in uint8 intensity units. + + ``contour_count`` : int + Number of external contours detected from ``phase_mask_binary``. + + Notes + ----- + The function applies standard non-periodic image connectivity. Hence, a + physical region crossing a periodic simulation boundary may be assigned + separate labels at opposite image edges. Periodic label merging should be + implemented before using region labels for periodic material identification or + periodic adaptive-grid decisions. + + The returned ``phase_mask_label`` describes disconnected geometric regions, + not material classes. In a future extension, regions with similar mean + grayscale values can be grouped into a separate ``material_label`` field, + while retaining ``phase_mask_label`` for region geometry. + """ + + _validate_odd_positive("blur_ksize", blur_ksize) + _validate_odd_positive("morph_kernel_size", morph_kernel_size) + + if blur_sigma < 0: + raise ValueError(f"blur_sigma must be non-negative, got {blur_sigma}.") + if open_iterations < 0: + raise ValueError( + f"open_iterations must be non-negative, got {open_iterations}." + ) + if close_iterations < 0: + raise ValueError( + f"close_iterations must be non-negative, got {close_iterations}." + ) + if connectivity not in (4, 8): + raise ValueError( + f"connectivity must be 4 or 8, got {connectivity}." + ) + + image_norm, image_u8 = _normalize_to_u8(data) + + image_blur = cv2.GaussianBlur( + image_u8, + (blur_ksize, blur_ksize), + blur_sigma, + ) + + otsu_threshold, mask_raw_u8 = cv2.threshold( + image_blur, + 0, + 255, + cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU, + ) + + kernel = np.ones( + (morph_kernel_size, morph_kernel_size), + dtype=np.uint8, + ) + + mask_open_u8 = cv2.morphologyEx( + mask_raw_u8, + cv2.MORPH_OPEN, + kernel, + iterations=open_iterations, + ) + + mask_clean_u8 = cv2.morphologyEx( + mask_open_u8, + cv2.MORPH_CLOSE, + kernel, + iterations=close_iterations, + ) + + # Binary phase indicator: values are always 0 or 1. + phase_mask_binary = (mask_clean_u8 > 0).astype(np.uint8) + + # Detect edges from the binary phase mask. + contours, _ = cv2.findContours( + (phase_mask_binary * 255).astype(np.uint8), + cv2.RETR_EXTERNAL, + cv2.CHAIN_APPROX_SIMPLE, + ) + + edge_mask = np.zeros_like(phase_mask_binary, dtype=np.uint8) + + cv2.drawContours( + edge_mask, + contours, + contourIdx=-1, + color=1, + thickness=1, + ) + + # Default behavior: label field is the same as the binary mask. + phase_mask_label = phase_mask_binary.copy() + + ''' + (1) connectedComponentsWithStats 對 binary foreground 做連通區域標記，得到每個 region 的唯一 ID。 + (2) 對每個 region，從原始 normalized grayscale image計算平均灰階。 + (3) 比較各 region 的平均灰階；平均值足夠相近的 region，給予相同 material_label。 + --> 會同時得到「每一塊幾何區域」與「每一種材料」兩種不同資料。connectedComponentsWithStats 可輸出 label map、每個區域統計與中心位置，而 label 0 是背景 + + [output]: + phase_mask_binary : 0/1，Otsu 分割結果 + phase_mask_label : 0/1/2/3/...，每個連通區域有唯一 ID + material_label : 0/1/2/...，平均灰階相近的區域共享材料 ID + region_mean_gray : 每個 region 的平均灰階 + region_to_material: region ID -> material ID 的對照表 + + def _validate_gray_tolerance(gray_tolerance: float) -> None: + """Validate the tolerance used to group region mean grayscale values.""" + if gray_tolerance < 0.0: + raise ValueError( + "gray_tolerance must be non-negative, " + f"got {gray_tolerance}." + ) + + + def _group_regions_by_mean_gray( + image_norm: np.ndarray, + phase_mask_label: np.ndarray, + gray_tolerance: float, + ) -> tuple[np.ndarray, dict[int, float], dict[int, int]]: + """ + Assign material IDs by grouping connected regions with similar mean gray. + + Parameters + ---------- + image_norm + Original image normalized to [0, 1]. + phase_mask_label + Connected-region label map. Label 0 is background. + gray_tolerance + Two regions are assigned to the same material when the absolute + difference between their representative mean gray values is not + greater than this tolerance. + + Returns + ------- + material_label + Integer material-ID field. Background remains 0. Foreground material + IDs start at 1. + + region_mean_gray + Dictionary mapping region ID to its mean normalized grayscale value. + + region_to_material + Dictionary mapping region ID to its assigned material ID. + """ + _validate_gray_tolerance(gray_tolerance) + + labels = np.asarray(phase_mask_label, dtype=np.int32) + material_label = np.zeros_like(labels, dtype=np.int32) + + region_ids = np.unique(labels) + region_ids = region_ids[region_ids != 0] + + region_mean_gray: dict[int, float] = {} + + for region_id in region_ids: + region_id = int(region_id) + region_pixels = image_norm[labels == region_id] + + if region_pixels.size == 0: + continue + + region_mean_gray[region_id] = float(region_pixels.mean()) + + # 每一組儲存： + # [代表平均灰階, 此組包含幾個 region] + material_groups: list[list[float]] = [] + + region_to_material: dict[int, int] = {} + + # 先依平均灰階排序，避免 label 的編號順序影響分類結果。 + sorted_region_ids = sorted( + region_mean_gray, + key=lambda region_id: region_mean_gray[region_id], + ) + + for region_id in sorted_region_ids: + region_gray = region_mean_gray[region_id] + material_id: int | None = None + + for group_index, group in enumerate(material_groups): + group_mean_gray = group[0] + group_region_count = group[1] + + if abs(region_gray - group_mean_gray) <= gray_tolerance: + material_id = group_index + 1 + + # 以目前 group 成員的平均值更新代表灰階。 + group[0] = ( + group_mean_gray * group_region_count + region_gray + ) / (group_region_count + 1) + + group[1] = group_region_count + 1 + break + + # 找不到相近 group，建立一個新材料類別。 + if material_id is None: + material_groups.append([region_gray, 1.0]) + material_id = len(material_groups) + + region_to_material[region_id] = material_id + material_label[labels == region_id] = material_id + + return material_label, region_mean_gray, region_to_material + ''' + # Optional: assign one integer label to each connected region. + if regions_label: + number_of_labels, phase_mask_label, phase_region_stats, phase_region_centroids = ( + cv2.connectedComponentsWithStats( + phase_mask_binary, + connectivity=connectivity, + ltype=cv2.CV_32S, + ) + ) + + results: dict[str, Any] = { + "edge_mask": edge_mask, + "phase_mask_binary": phase_mask_binary, + "phase_mask_label": phase_mask_label, + } + + # These data exist only when connected-component labeling was requested. + if regions_label: + results.update( + { + "number_of_phase_regions": int(number_of_labels - 1), + "phase_region_stats": phase_region_stats, + "phase_region_centroids": phase_region_centroids, + } + ) + + if return_intermediate: + results.update( + { + "image_norm": image_norm, + "image_u8": image_u8, + "image_blur": image_blur, + "mask_raw": (mask_raw_u8 > 0).astype(np.uint8), + "mask_open": (mask_open_u8 > 0).astype(np.uint8), + "mask_clean": phase_mask_binary, + "otsu_threshold": float(otsu_threshold), + "contour_count": int(len(contours)), + } + ) + + return results