From 66edc029d21d301ecda600e64f787a5b5386dc3f Mon Sep 17 00:00:00 2001 From: MartinLadecky Date: Tue, 30 Jun 2026 16:18:43 +0200 Subject: [PATCH 1/3] Adding examples of homogenization on a adapted grid. --- ...ation_conductivity_transformed_grid_Jia.py | 212 +++ ...ization_elasticity_transformed_grid_Jia.py | 233 +++ .../exp_analytical_grid_adaptation.py | 166 +++ .../exp_circle_muGrid_real_field.py | 340 +++++ muFFTTO/analytical_grid_adaptation.py | 479 +++++++ muFFTTO/grid_adaptation_methods.py | 1245 +++++++++++++++++ 6 files changed, 2675 insertions(+) create mode 100644 experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia.py create mode 100644 experiments/grid_adaptation/exp_2D_homogenization_elasticity_transformed_grid_Jia.py create mode 100644 experiments/grid_adaptation/exp_analytical_grid_adaptation.py create mode 100644 experiments/grid_adaptation/exp_circle_muGrid_real_field.py create mode 100644 muFFTTO/analytical_grid_adaptation.py create mode 100644 muFFTTO/grid_adaptation_methods.py 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_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/analytical_grid_adaptation.py b/muFFTTO/analytical_grid_adaptation.py new file mode 100644 index 0000000..06a823b --- /dev/null +++ b/muFFTTO/analytical_grid_adaptation.py @@ -0,0 +1,479 @@ +''' +Storage of adapted grid data into muGrid real_field containers +through the primary accessor `field.p`. + +Reference: +Zecevic, M., Lebensohn, R. A., & Capolungo, L. (2026). +Achieving geometric accuracy in FFT-based micromechanical models +using conformal grid. Mechanics of Materials, 212, 105512. +https://doi.org/10.1016/j.mechmat.2025.105512 + +''' + +import numpy as np +from collections import deque +import muGrid + +def make_grid_nodes(N: int = 64, L: float = 25.0) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Function that constructs a regular Cartesian node grid. + + The grid has (N+1) x (N+1) nodes that span a square domain [-L, L] x [-L, L]. + Node coordinates are stored as an array of shape [nx, ny, xy]. + + Parameters + -------- + N : int + Number of cells along one spatial direction. The number of nodes is N + 1. + L : float + Half-size of the square computational domain measured in the same units as x and y. + + Returns + ------- + [i,x,y] (x,y),nx,ny + P : numpy ndarray + Array of nodal coordinates with shape [xy, nx, ny]. + - nx = N is the number of nodes in x-direction. + - ny = N is the number of nodes in y-direction. + - the last index 'xy' holds (x, y) coordinates at each node. + x : numpy ndarray + One-dimensional array of x-coordinates with shape [nx]. + y : numpy ndarray + One-dimensional array of y-coordinates with shape [ny]. + """ + + x = np.linspace(-L, L, N+1) + y = np.linspace(-L, L, N+1) + X, Y = np.meshgrid(x, y, indexing="ij") + P = np.stack([X, Y], axis=-1) + + ''' + # (0,L) + x = np.linspace(0, L, N, endpoint=False ) + y = np.linspace(0, L, N, endpoint=False ) + X, Y = np.meshgrid(x, y, indexing="ij") + P = np.stack([X, Y]) + ''' + return P, x, y + + +# specifically for circle inclusion, but can be adapted for other shapes by changing the labeling logic +def cell_labels( + P: np.ndarray, + center: tuple[float, float] = (0.0, 0.0), + R: float = 20.0, +) -> np.ndarray: + """ + Function that labels grid cells as inside or outside a circular inclusion. + + The label is decided from the position of cell centers with respect to a circle. + A value of one means inside the circle, zero means outside. + + Parameters + ---------- + P : numpy ndarray + Nodal coordinates with shape [nx, ny, xy] as returned by make_grid_nodes. + - nx, ny are nodal counts in x- and y-direction. + - xy index stores (x, y) coordinates per node. + center : tuple of float + Coordinates of the circle center (cx, cy). + R : float + Radius of the inclusion circle. + + Returns + ------- + inside : numpy ndarray + Integer array of cell labels with shape [cell_x, cell_y]. + - inside[i, j] = 1 if the cell center lies inside the circle. + - inside[i, j] = 0 otherwise. + """ + cx, cy = center + # cell centers from average of four surrounding nodes + Pc = 0.25 * (P[:-1, :-1] + P[1:, :-1] + P[:-1, 1:] + P[1:, 1:]) + dx = Pc[..., 0] - cx + dy = Pc[..., 1] - cy + inside_mask = (dx * dx + dy * dy) <= R * R + return inside_mask.astype(np.int32) + + +def interface_node_mask(cell_inside: np.ndarray) -> np.ndarray: + """ + Function that detects interface nodes from neighboring cell labels. + + A node is classified as interface node if the adjacent cells contain both inside and outside labels. + + Parameters + ---------- + cell_inside : numpy ndarray + Cell labels with shape [cell_x, cell_y] created by cell_labels. + - value 1 corresponds to cells inside the circle. + - value 0 corresponds to cells outside the circle. + + Returns + ------- + interface_mask : numpy ndarray + Boolean nodal mask with shape [nx, ny]. + - True indicates that the node lies on the interface between phases. + - False indicates a node far from the interface. + """ + N = cell_inside.shape[0] + mask = np.zeros((N + 1, N + 1), dtype=bool) + + for i in range(N + 1): + for j in range(N + 1): + values = [] + if i > 0 and j > 0: + values.append(cell_inside[i - 1, j - 1]) + if i < N and j > 0: + values.append(cell_inside[i, j - 1]) + if i > 0 and j < N: + values.append(cell_inside[i - 1, j]) + if i < N and j < N: + values.append(cell_inside[i, j]) + if not values: + continue + vals = np.asarray(values, dtype=int) + if vals.size >= 2 and vals.min() != vals.max(): + mask[i, j] = True + + return mask + + +def outer_boundary_mask(N: int) -> np.ndarray: + """ + Function that marks the outer square boundary nodes as fixed. + + The boundary corresponds to the four edges of the structured grid. + The mask can be used to impose Dirichlet-type constraints in relaxation. + + Parameters + ---------- + N : int + Number of cells along one spatial direction. The number of nodes is N + 1. + + Returns + ------- + outer_mask : numpy ndarray + Boolean nodal mask with shape [nx, ny]. + - True indicates a node on the outer boundary of the domain. + - False indicates an interior node. + """ + m = np.zeros((N + 1, N + 1), dtype=bool) + m[0, :] = True + m[N, :] = True + m[:, 0] = True + m[:, N] = True + return m + +#specifically for circle inclusion, but can be adapted for other shapes by changing the projection logic +def project_points_to_circle( + Ppts: np.ndarray, + center: tuple[float, float] = (0.0, 0.0), + R: float = 20.0, +) -> np.ndarray: + """ + Function that projects given points onto a target circle. + + Each point is shifted radially such that its distance to the center + equals the prescribed radius R. + + Parameters + ---------- + Ppts : numpy ndarray + Array of point coordinates with shape [k, xy]. + - k is the number of points to project. + - xy index stores (x, y) coordinates per point. + center : tuple of float + Coordinates (cx, cy) of the circle center. + R : float + Target radius of the circle. + + Returns + ------- + Pproj : numpy ndarray + Array of projected point coordinates with shape [k, xy]. + - all returned points satisfy the distance ||Pproj - center|| = R within numerical tolerance. + """ + c = np.asarray(center, dtype=float) + V = Ppts - c[None, :] + r = np.linalg.norm(V, axis=1, keepdims=True) + r = np.maximum(r, 1e-12) + return c[None, :] + (R / r) * V + + +def manhattan_distance_to_interface(interface_mask: np.ndarray) -> np.ndarray: + """ + Function that computes Manhattan distance to the nearest interface node. + + The distance is measured on the grid graph using unit steps in the + four von Neumann directions (up, down, left, right). + + Parameters + ---------- + interface_mask : numpy ndarray + Boolean nodal mask with shape [nx, ny]. + - True marks nodes that belong to the interface. + - False marks other nodes. + + Returns + ------- + dist : numpy ndarray + Floating-point distance field with shape [nx, ny]. + - dist[i, j] = 0 for interface nodes. + - dist[i, j] is the minimum number of grid steps from node (i, j) to the interface. + """ + N = interface_mask.shape[0] - 1 + dist = np.full((N + 1, N + 1), np.inf, dtype=float) + queue = deque() + + for i in range(N + 1): + for j in range(N + 1): + if interface_mask[i, j]: + dist[i, j] = 0.0 + queue.append((i, j)) + + neighbors = [(1, 0), (-1, 0), (0, 1), (0, -1)] + while queue: + i, j = queue.popleft() + for di, dj in neighbors: + ii, jj = i + di, j + dj + if 0 <= ii <= N and 0 <= jj <= N: + if dist[ii, jj] > dist[i, j] + 1.0: + dist[ii, jj] = dist[i, j] + 1.0 + queue.append((ii, jj)) + + return dist + + +def stiffness_from_distance( + dist: np.ndarray, + k0: float = 1.0, + b: float = 1.0, + a: float = 1.0, + kmin: float = 1e-3, +) -> np.ndarray: + """ + Function that computes nodal stiffness from distance to the interface. + + The stiffness follows a decay law k(d) = k0 / (d + a)^b and is bounded from below + by the stiffness floor kmin. + + Parameters + ---------- + dist : numpy ndarray + Manhattan distance field with shape [nx, ny]. + - dist[i, j] is the distance of node (i, j) to the interface. + k0 : float + Baseline stiffness scale factor. + b : float + Exponent that controls how fast stiffness decays away from the interface. + a : float + Positive offset that regularizes the stiffness at zero distance. + kmin : float + Minimum stiffness floor applied elementwise. + + Returns + ------- + k_node : numpy ndarray + Nodal stiffness field with shape [nx, ny]. + - k_node[i, j] is the effective stiffness assigned to node (i, j). + - values are guaranteed to be >= kmin. + """ + k = k0 / np.power(dist + a, b) + return np.maximum(k, kmin) + + +def spring_relax_weighted( + P: np.ndarray, + fixed_mask: np.ndarray, + k_node: np.ndarray, + iters: int = 600, + omega: float = 1.0, +) -> np.ndarray: + """ + Function that performs weighted spring-based relaxation of nodal positions. + + The update is a Gauss–Seidel type iteration that moves each free node + toward the weighted average of its four neighbors, where edge weights + are computed from nodal stiffness values. + + Parameters + ---------- + P : numpy ndarray + Nodal coordinates before relaxation with shape [nx, ny, xy]. + fixed_mask : numpy ndarray + Boolean nodal mask with shape [nx, ny] that marks fixed nodes. + - True indicates that the node is held fixed during the update. + k_node : numpy ndarray + Nodal stiffness field with shape [nx, ny]. + - stiffness values control how strongly a node couples to its neighbors. + iters : int + Number of Gauss–Seidel sweeps over all interior nodes. + omega : float + Relaxation parameter. Values smaller than one yield under-relaxation. + + Returns + ------- + P_relaxed : numpy ndarray + Nodal coordinates after relaxation with shape [nx, ny, xy]. + - fixed nodes remain unchanged. + - free nodes are updated to a spring-equilibrium configuration. + """ + P_new = P.copy() + N = P.shape[0] - 1 + + for _ in range(iters): + for i in range(1, N): + for j in range(1, N): + if fixed_mask[i, j]: + continue + + neighbors = [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)] + kij = k_node[i, j] + w = np.empty(4, dtype=float) + pts = np.empty((4, 2), dtype=float) + + for idx, (ii, jj) in enumerate(neighbors): + w[idx] = 0.5 * (kij + k_node[ii, jj]) + pts[idx] = P_new[ii, jj] + + target = (w[:, None] * pts).sum(axis=0) / (w.sum() + 1e-15) + P_new[i, j] = (1.0 - omega) * P_new[i, j] + omega * target + + return P_new + + +def adapt_grid_to_circle( + N: int = 64, + L: float = 40.0, + center: tuple[float, float] = (0.0, 0.0), + R: float = 20.0, + iters: int = 600, + omega: float = 0.8, + b: float = 1.0, + k0: float = 0.25, + kmin: float = 0.005, +) -> dict: + """ + Function that adapts a structured grid to a circular inclusion. + + The algorithm performs cell labeling, interface detection, projection + of interface nodes onto the circle and weighted spring relaxation of + interior nodes. + + Parameters + ---------- + N : int + Number of cells along one spatial direction. + L : float + Half-size of the square computational domain. + center : tuple of float + Coordinates (cx, cy) of the circle center. + R : float + Radius of the inclusion circle. + iters : int + Number of relaxation sweeps in spring_relax_weighted. + omega : float + Relaxation parameter for spring_relax_weighted. + b : float + Exponent of the distance-based stiffness decay law. + k0 : float + Baseline stiffness parameter. + kmin : float + Minimum stiffness floor. + + Returns + ------- + result : dict + Dictionary that collects all relevant fields: + - 'P0' : initial nodal coordinates [nx, ny, xy] + - 'P' : adapted nodal coordinates [nx, ny, xy] + - 'inside' : cell labels [cell_x, cell_y] + - 'interface' : interface nodal mask [nx, ny] + - 'outer' : outer boundary mask [nx, ny] + - 'fixed' : fixed nodal mask [nx, ny] + - 'dist' : Manhattan distance field [nx, ny] + - 'k_node' : nodal stiffness values [nx, ny] + - 'params' : dictionary of scalar parameters used in the run. + """ + P0, _, _ = make_grid_nodes(N=N, L=L) + inside = cell_labels(P0, center=center, R=R) + interface = interface_node_mask(inside) + outer = outer_boundary_mask(N) + fixed = outer | interface + + P1 = P0.copy() + idx = np.argwhere(interface) + if idx.size: + P_if = P1[idx[:, 0], idx[:, 1], :] + P1[idx[:, 0], idx[:, 1], :] = project_points_to_circle( + P_if, + center=center, + R=R, + ) + + dist = manhattan_distance_to_interface(interface) + k_node = stiffness_from_distance(dist, k0=k0, b=b, a=1.0, kmin=kmin) + P2 = spring_relax_weighted(P1, fixed_mask=fixed, k_node=k_node, iters=iters, omega=omega) + + return { + "P0": P0, + "P": P2, + "inside": inside, + "interface": interface, + "outer": outer, + "fixed": fixed, + "dist": dist, + "k_node": k_node, + "params": { + "N": N, + "L": L, + "center": center, + "R": R, + "iters": iters, + "omega": omega, + "b": b, + "k0": k0, + "kmin": kmin, + }, + } + + + +''' +from muGrid import Communicator, CartesianDecomposition +import numpy as np + +if __name__ == "__main__": + coor,x,z = make_grid_nodes(N=4, L=1.0) + print(coor[0,0]) + print(coor[1,0]) + + #print(muGrid.__file__) + #print([x for x in dir(muGrid) if not x.startswith("_")]) + + nx = 8 + ny = 8 + g = 1 + print("before communicator") + comm = Communicator() + print("after communicator") + print("before decomp") + + #help(muGrid.CartesianDecomposition)vu0 + decomp = muGrid.CartesianDecomposition( + communicator=comm, + nb_domain_grid_pts=(nx, ny), + nb_subdivisions=(1, 1), + nb_ghosts_left=(g, g), + nb_ghosts_right=(g, g), + ) + print("after decomp") + + field = decomp.real_field("displacement", components=(3,)) + print(comm) + print(field) + print(type(field)) + print([x for x in dir(field) if not x.startswith("_")]) + +''' \ No newline at end of file diff --git a/muFFTTO/grid_adaptation_methods.py b/muFFTTO/grid_adaptation_methods.py new file mode 100644 index 0000000..4f5829f --- /dev/null +++ b/muFFTTO/grid_adaptation_methods.py @@ -0,0 +1,1245 @@ +""" +Storage into muGrid +real_field containers through the primary accessor `field.p`. + +Reference: +Zecevic, M., Lebensohn, R. A., & Capolungo, L. (2026). +Achieving geometric accuracy in FFT-based micromechanical models +using conformal grid. Mechanics of Materials, 212, 105512. +https://doi.org/10.1016/j.mechmat.2025.105512 +""" + +from __future__ import annotations +from collections import deque +import numpy as np +import muGrid + + +def make_grid_nodes(N: int = 64, L: float = 25.0) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Function that constructs a regular periodic-style Cartesian node grid. + + The grid stores N x N nodal points on the square domain [0, L) x [0, L). + The right and top boundary nodes are not stored explicitly, which is consistent + with a periodic grid representation. Node coordinates are stored in an array + of shape [xy, nx, ny]. + + Parameters + ---------- + N : int + Number of stored grid points per spatial direction. + L : float + Side length of the square computational domain measured in the same units + as x and y. + + Returns + ------- + P : numpy.ndarray + Array of nodal coordinates with shape [xy, nx, ny]. + - xy = 2 stores the x- and y-coordinates. + - nx = N is the number of stored nodes in x-direction. + - ny = N is the number of stored nodes in y-direction. + x : numpy.ndarray + One-dimensional array of x-coordinates with shape [nx]. + y : numpy.ndarray + One-dimensional array of y-coordinates with shape [ny]. + """ + x = np.linspace(0.0, L, N, endpoint=False) + y = np.linspace(0.0, L, N, endpoint=False) + X, Y = np.meshgrid(x, y, indexing="ij") + P = np.stack([X, Y], axis=0) + return P, x, y + + +def cell_labels( + P: np.ndarray, + center: tuple[float, float] = (0.0, 0.0), + R: float = 20.0, +) -> np.ndarray: + """ + Function that labels cells as inside or outside a circular inclusion. + + Cell centers are computed from the four corner nodes of each periodic cell. + A cell is marked as inside if its center lies inside or on the target circle. + + Parameters + ---------- + P : numpy.ndarray + Array of nodal coordinates with shape [xy, nx, ny]. + center : tuple[float, float] + Coordinates of the circle center given as (x_center, y_center). + R : float + Radius of the circular inclusion. + + Returns + ------- + inside : numpy.ndarray + Integer array with shape [nx, ny]. + A value of 1 indicates that the periodic cell center lies inside the circle, + and a value of 0 indicates that it lies outside. + """ + cx, cy = center + + P_ip1 = np.roll(P, shift=-1, axis=1) + P_jp1 = np.roll(P, shift=-1, axis=2) + P_ip1_jp1 = np.roll(P_ip1, shift=-1, axis=2) + + P_ip1 = P_ip1.copy() + P_jp1 = P_jp1.copy() + P_ip1_jp1 = P_ip1_jp1.copy() + + dx_grid = P[0, 1, 0] - P[0, 0, 0] + dy_grid = P[1, 0, 1] - P[1, 0, 0] + + P_ip1[0, -1, :] += dx_grid * P.shape[1] + P_ip1_jp1[0, -1, :] += dx_grid * P.shape[1] + P_jp1[1, :, -1] += dy_grid * P.shape[2] + P_ip1_jp1[1, :, -1] += dy_grid * P.shape[2] + + Pc = 0.25 * (P + P_ip1 + P_jp1 + P_ip1_jp1) + dx = Pc[0] - cx + dy = Pc[1] - cy + inside_mask = (dx * dx + dy * dy) <= R * R + return inside_mask.astype(np.int32) + + +def interface_node_mask(cell_inside: np.ndarray) -> np.ndarray: + """ + Function that detects interface nodes on a periodic stored grid. + + A stored node is marked as an interface node if the surrounding periodic cells + do not all have the same inside/outside label. In other words, the node lies + on the discrete boundary between two material regions. + + Parameters + ---------- + cell_inside : numpy.ndarray + Integer array of cell labels with shape [nx, ny]. + Cells with value 1 are inside the inclusion and cells with value 0 are outside. + + Returns + ------- + mask : numpy.ndarray + Boolean array with shape [nx, ny]. + True indicates that the node is an interface node, and False otherwise. + """ + N = cell_inside.shape[0] + mask = np.zeros((N, N), dtype=bool) + + for i in range(N): + for j in range(N): + vals = np.array([ + cell_inside[(i - 1) % N, (j - 1) % N], + cell_inside[i % N, (j - 1) % N], + cell_inside[(i - 1) % N, j % N], + cell_inside[i % N, j % N], + ], dtype=int) + if vals.min() != vals.max(): + mask[i, j] = True + + return mask + + +def project_points_to_circle( + Ppts: np.ndarray, + center: tuple[float, float] = (0.0, 0.0), + R: float = 20.0, +) -> np.ndarray: + """ + Function that projects points onto a target circle. + + Each input point is moved along the radial direction so that its distance + from the circle center becomes exactly equal to the prescribed radius. + The function accepts point arrays in either [xy, k] or [k, xy] format. + + Parameters + ---------- + Ppts : numpy.ndarray + Array of point coordinates with shape [xy, k] or [k, xy]. + center : tuple[float, float] + Coordinates of the circle center given as (x_center, y_center). + R : float + Radius of the target circle. + + Returns + ------- + proj : numpy.ndarray + Array of projected point coordinates with the same shape convention + as the input array `Ppts`. + + Raises + ------ + ValueError + If `Ppts` is not a two-dimensional array with shape [xy, k] or [k, xy]. + """ + c = np.asarray(center, dtype=float) + + transposed = False + if Ppts.ndim != 2: + raise ValueError("Ppts must be a 2D array") + if Ppts.shape[0] == 2: + pts = Ppts.T + transposed = True + elif Ppts.shape[1] == 2: + pts = Ppts + else: + raise ValueError("Ppts must have shape [xy, k] or [k, xy]") + + V = pts - c[None, :] + r = np.linalg.norm(V, axis=1, keepdims=True) + r = np.maximum(r, 1e-12) + proj = c[None, :] + (R / r) * V + return proj.T if transposed else proj + + +def manhattan_distance_to_interface(interface_mask: np.ndarray) -> np.ndarray: + """ + Function that computes the periodic Manhattan distance to the nearest interface node. + + The distance is measured in the discrete grid sense using nearest-neighbor + connectivity in the x- and y-directions. Periodic wrapping is applied at the + domain boundaries. + + Parameters + ---------- + interface_mask : numpy.ndarray + Boolean array with shape [nx, ny]. + True marks interface nodes and False marks non-interface nodes. + + Returns + ------- + dist : numpy.ndarray + Floating-point array with shape [nx, ny]. + Each entry contains the periodic Manhattan distance from that node to + the nearest interface node. + """ + N = interface_mask.shape[0] + dist = np.full((N, N), np.inf, dtype=float) + queue = deque() + + for i in range(N): + for j in range(N): + if interface_mask[i, j]: + dist[i, j] = 0.0 + queue.append((i, j)) + + neighbors = [(1, 0), (-1, 0), (0, 1), (0, -1)] + while queue: + i, j = queue.popleft() + for di, dj in neighbors: + ii = (i + di) % N + jj = (j + dj) % N + if dist[ii, jj] > dist[i, j] + 1.0: + dist[ii, jj] = dist[i, j] + 1.0 + queue.append((ii, jj)) + + return dist + + +def stiffness_from_distance( + dist: np.ndarray, + k0: float = 1.0, + b: float = 1.0, + a: float = 1.0, + kmin: float = 1e-3, +) -> np.ndarray: + """ + Function that computes a nodal stiffness field from the distance to the interface. + + The stiffness is defined by a distance-dependent power law and is bounded + from below by a prescribed minimum stiffness value. + + Parameters + ---------- + dist : numpy.ndarray + Array of nodal distances to the nearest interface with shape [nx, ny]. + k0 : float + Reference stiffness factor. + b : float + Exponent controlling how fast the stiffness decays with distance. + a : float + Positive shift added to the distance to avoid division by zero and to + control the stiffness near the interface. + kmin : float + Minimum admissible stiffness value. + + Returns + ------- + k : numpy.ndarray + Floating-point stiffness array with shape [nx, ny]. + """ + k = k0 / np.power(dist + a, b) + return np.maximum(k, kmin) + + +def spring_relax_weighted( + P: np.ndarray, + fixed_mask: np.ndarray, + k_node: np.ndarray, + iters: int = 600, + omega: float = 1.0, +) -> np.ndarray: + """ + Function that performs weighted spring relaxation on a periodic stored grid. + + Each non-fixed node is iteratively moved toward a weighted average of its + four nearest neighbors. The nodal stiffness field controls the local strength + of the spring interaction, and periodic wrapping is applied at the domain boundaries. + + Parameters + ---------- + P : numpy.ndarray + Array of nodal coordinates with shape [xy, nx, ny]. + fixed_mask : numpy.ndarray + Boolean array with shape [nx, ny]. + True marks fixed nodes that remain unchanged during relaxation. + k_node : numpy.ndarray + Floating-point nodal stiffness array with shape [nx, ny]. + iters : int + Number of relaxation iterations. + omega : float + Relaxation parameter. Values between 0 and 1 correspond to under-relaxation, + and omega = 1 applies the full update at each iteration. + + Returns + ------- + P_new : numpy.ndarray + Relaxed nodal coordinates with shape [xy, nx, ny]. + """ + P_new = P.copy() + N = P.shape[1] + + dx_grid = P[0, 1, 0] - P[0, 0, 0] + dy_grid = P[1, 0, 1] - P[1, 0, 0] + + for _ in range(iters): + for i in range(N): + for j in range(N): + if fixed_mask[i, j]: + continue + + neighbors = [ + ((i - 1) % N, j), + ((i + 1) % N, j), + (i, (j - 1) % N), + (i, (j + 1) % N), + ] + kij = k_node[i, j] + w = np.empty(4, dtype=float) + pts = np.empty((4, 2), dtype=float) + + xij = P_new[0, i, j] + yij = P_new[1, i, j] + + for idx, (ii, jj) in enumerate(neighbors): + xnb = P_new[0, ii, jj] + ynb = P_new[1, ii, jj] + + if ii == 0 and i == N - 1: + xnb += dx_grid * N + elif ii == N - 1 and i == 0: + xnb -= dx_grid * N + + if jj == 0 and j == N - 1: + ynb += dy_grid * N + elif jj == N - 1 and j == 0: + ynb -= dy_grid * N + + w[idx] = 0.5 * (kij + k_node[ii, jj]) + pts[idx] = [xnb, ynb] + + target = (w[:, None] * pts).sum(axis=0) / (w.sum() + 1e-15) + P_new[:, i, j] = (1.0 - omega) * np.array([xij, yij]) + omega * target + + return P_new + + +def adapt_grid_to_circle( + N: int = 64, + L: float = 40.0, + center: tuple[float, float] = (20.0, 20.0), + R: float = 10.0, + iters: int = 600, + omega: float = 0.8, + b: float = 1.0, + k0: float = 0.25, + kmin: float = 0.005, +) -> dict: + """ + Function that adapts a periodic Cartesian grid to a circular inclusion. + + The procedure consists of four main steps: + 1. Construct the initial periodic grid. + 2. Label cells and detect interface nodes. + 3. Project interface nodes onto the target circle. + 4. Relax the remaining nodes with a weighted spring model. + + Parameters + ---------- + N : int + Number of stored grid points per spatial direction. + L : float + Side length of the square computational domain. + center : tuple[float, float] + Coordinates of the circle center given as (x_center, y_center). + R : float + Radius of the circular inclusion. + iters : int + Number of spring-relaxation iterations. + omega : float + Relaxation parameter for the spring-relaxation update. + b : float + Exponent of the stiffness-distance law. + k0 : float + Reference stiffness factor. + kmin : float + Minimum admissible stiffness value. + + Returns + ------- + result : dict + Dictionary containing the adapted grid and related fields: + - "P0": initial nodal coordinates, shape [xy, nx, ny] + - "P": adapted nodal coordinates, shape [xy, nx, ny] + - "inside": cell labels, shape [nx, ny] + - "interface": interface-node mask, shape [nx, ny] + - "fixed": fixed-node mask, shape [nx, ny] + - "dist": distance-to-interface field, shape [nx, ny] + - "k_node": nodal stiffness field, shape [nx, ny] + - "params": dictionary of input parameters + """ + P0, _, _ = make_grid_nodes(N=N, L=L) + inside = cell_labels(P0, center=center, R=R) + interface = interface_node_mask(inside) + fixed = interface.copy() + + P1 = P0.copy() + idx = np.argwhere(interface) + if idx.size: + pts = P1[:, idx[:, 0], idx[:, 1]] + P1[:, idx[:, 0], idx[:, 1]] = project_points_to_circle( + pts, + center=center, + R=R, + ) + + dist = manhattan_distance_to_interface(interface) + k_node = stiffness_from_distance(dist, k0=k0, b=b, a=1.0, kmin=kmin) + P2 = spring_relax_weighted(P1, fixed_mask=fixed, k_node=k_node, iters=iters, omega=omega) + + return { + "P0": P0, + "P": P2, + "inside": inside, + "interface": interface, + "fixed": fixed, + "dist": dist, + "k_node": k_node, + "params": { + "N": N, + "L": L, + "center": center, + "R": R, + "iters": iters, + "omega": omega, + "b": b, + "k0": k0, + "kmin": kmin, + }, + } + + +def adapt_grid_to_circle_EXAMPLE_( + nb_grid_points: tuple[int, int] = (64, 64), + domain_size: tuple[float, float] = (1, 1), + center: tuple[float, float] = (0.5, 0.5), + radius: float = 0.2, + reference_grid_points_coords: np.ndarray = None, + iters: int = 600, + omega: float = 0.8, + b: float = 1.0, + k0: float = 0.25, + kmin: float = 0.005, +) -> dict: + """ + Function that adapts a periodic Cartesian grid to a circular inclusion. + + The procedure consists of four main steps: + 1. Construct the initial periodic grid. + 2. Label cells and detect interface nodes. + 3. Project interface nodes onto the target circle. + 4. Relax the remaining nodes with a weighted spring model. + + Parameters + ---------- + nb_grid_points : tuple[int, int] + Number of stored grid points per spatial direction. + domain_size : tuple[float, float] + Side lengths of the rectangular computational domain. + center : tuple[float, float] + Coordinates of the circle center given as (x_center, y_center). + radius : float + Radius of the circular inclusion. + iters : int + Number of spring-relaxation iterations. + omega : float + Relaxation parameter for the spring-relaxation update. + b : float + Exponent of the stiffness-distance law. + k0 : float + Reference stiffness factor. + kmin : float + Minimum admissible stiffness value. + + Returns + ------- + result : dict + Dictionary containing the adapted grid and related fields: + - "P0": initial nodal coordinates, shape [xy, nx, ny] + - "P": adapted nodal coordinates, shape [xy, nx, ny] + - "inside": cell labels, shape [nx, ny] + - "interface": interface-node mask, shape [nx, ny] + - "fixed": fixed-node mask, shape [nx, ny] + - "dist": distance-to-interface field, shape [nx, ny] + - "k_node": nodal stiffness field, shape [nx, ny] + - "params": dictionary of input parameters + """ + P0 = reference_grid_points_coords # this is just aliasing, not copy + #P0, _, _ = make_grid_nodes(N=nb_grid_points[0], L=domain_size[0]) + inside = cell_labels(P0, center=center, R=radius) + interface = interface_node_mask(inside) + fixed = interface.copy() + + P1 = P0.copy() + idx = np.argwhere(interface) + if idx.size: + pts = P1[:, idx[:, 0], idx[:, 1]] + P1[:, idx[:, 0], idx[:, 1]] = project_points_to_circle( + pts, + center=center, + R=radius, + ) + # TODO[Jia]: what is this + dist = manhattan_distance_to_interface(interface) + # TODO[Jia]: Comment what is this + k_node = stiffness_from_distance(dist, k0=k0, b=b, a=1.0, kmin=kmin) + P1 = spring_relax_weighted(P1, fixed_mask=fixed, k_node=k_node, iters=iters, omega=omega) + + return { + "coords_of_displaced_nodes": P1, + "inside": inside, + "interface": interface, + "fixed": fixed, + "dist": dist, + "k_node": k_node, + "params": { + "nb_grid_points": nb_grid_points, + "domain_size": domain_size, + "center": center, + "radius": radius, + "iters": iters, + "omega": omega, + "b": b, + "k0": k0, + "kmin": kmin, + }, + } + + +def make_mugrid_decomposition( + nx: int, + ny: int, + ghosts: int = 1, +): + """ + Function that creates a simple single-process muGrid Cartesian decomposition. + + The decomposition is defined for a structured grid with a prescribed number + of domain points and a prescribed number of ghost points on each side. + + Parameters + ---------- + nx : int + Number of domain grid points in x-direction. + ny : int + Number of domain grid points in y-direction. + ghosts : int + Number of ghost points added on the left and right sides of each direction. + + Returns + ------- + comm : muGrid.Communicator + muGrid communicator object. + decomp : muGrid.CartesianDecomposition + Cartesian decomposition associated with the communicator and grid layout. + """ + comm = muGrid.Communicator() + decomp = muGrid.CartesianDecomposition( + communicator=comm, + nb_domain_grid_pts=(nx, ny), + nb_subdivisions=(1, 1), + nb_ghosts_left=(ghosts, ghosts), + nb_ghosts_right=(ghosts, ghosts), + ) + return comm, decomp + + +def _copy_array_into_field_p(field, arr: np.ndarray) -> None: + """ + Function that copies a NumPy array into a muGrid field through the primary accessor `p`. + + This function assumes that the muGrid field uses the shape convention + [components, nx, ny] for vector and scalar fields stored through `field.p`. + + Parameters + ---------- + field : muGrid.Field + muGrid field object whose primary accessor `p` will be written. + arr : numpy.ndarray + Array to be copied into the field. + Expected shapes are: + - [ncomp, nx, ny] for vector-valued fields + - [1, nx, ny] for scalar-valued fields stored with one component + + Returns + ------- + None + + Raises + ------ + ValueError + If the shape of `arr` does not match the shape of `field.p`. + """ + target = np.asarray(field.p) + if target.shape != arr.shape: + raise ValueError( + f"Shape mismatch for field.p: target shape = {target.shape}, arr shape = {arr.shape}" + ) + target[...] = arr + + +def make_numpy_field_bundle(result: dict, store_positions: bool = True) -> dict: + """ + Function that builds a dictionary of useful NumPy fields from the grid-adaptation result. + + The returned bundle contains deformation-related scalar and vector fields + derived from the output of `adapt_grid_to_circle`. The displacement field is + computed as the difference between the adapted and initial nodal coordinates. + + Parameters + ---------- + result : dict + Dictionary returned by `adapt_grid_to_circle`. + store_positions : bool + If True, the initial and adapted nodal coordinates `P0` and `P` are also stored + in the returned bundle. If False, only derived fields are stored. + + Returns + ------- + bundle : dict + Dictionary containing NumPy arrays: + - "displacement": shape [2, nx, ny] + - "inside": shape [1, nx, ny] + - "interface": shape [1, nx, ny] + - "fixed": shape [1, nx, ny] + - "dist": shape [1, nx, ny] + - "k_node": shape [1, nx, ny] + and optionally + - "P0": shape [2, nx, ny] + - "P": shape [2, nx, ny] + """ + P0 = result["P0"] + P = result["P"] + inside = result["inside"] + interface = result["interface"] + fixed = result["fixed"] + dist = result["dist"] + k_node = result["k_node"] + + displacement = P - P0 + + bundle = { + "displacement": displacement.astype(np.float64), + "inside": inside[None, ...].astype(np.float64), + "interface": interface[None, ...].astype(np.float64), + "fixed": fixed[None, ...].astype(np.float64), + "dist": dist[None, ...].astype(np.float64), + "k_node": k_node[None, ...].astype(np.float64), + } + + if store_positions: + bundle["P0"] = P0.astype(np.float64) + bundle["P"] = P.astype(np.float64) + + return bundle + + +def pack_adapted_grid_to_fields( + result: dict, + ghosts: int = 1, + store_positions: bool = True, + verbose: bool = True, +) -> dict: + """ + Function that stores adapted-grid data into muGrid real_field containers. + + The function first builds a NumPy field bundle from the adaptation result, + then creates one muGrid real_field per stored quantity, and finally copies + each NumPy array into the corresponding field through `field.p`. + + Parameters + ---------- + result : dict + Dictionary returned by `adapt_grid_to_circle`. + ghosts : int + Number of ghost points used in the muGrid decomposition. + store_positions : bool + If True, the initial and adapted nodal coordinates `P0` and `P` are also stored. + verbose : bool + If True, a short write-status message is printed for each stored field. + + Returns + ------- + bundle : dict + Dictionary containing: + - "comm": muGrid communicator + - "decomp": muGrid Cartesian decomposition + - "fields": dictionary of muGrid real_field objects + - "numpy": dictionary of NumPy arrays used for storage + - "report": dictionary summarizing shapes and storage status + """ + nx, ny = result["P"].shape[1], result["P"].shape[2] + comm, decomp = make_mugrid_decomposition(nx=nx, ny=ny, ghosts=ghosts) + + numpy_fields = make_numpy_field_bundle(result, store_positions=store_positions) + mugrid_fields = {} + report = {} + + for name, arr in numpy_fields.items(): + ncomp = arr.shape[0] if arr.ndim == 3 else 1 + field = decomp.real_field(name, components=(ncomp,)) + mugrid_fields[name] = field + + _copy_array_into_field_p(field, arr) + + report[name] = { + "status": "ok", + "shape": arr.shape, + "p_shape": np.asarray(field.p).shape, + "pg_shape": np.asarray(field.pg).shape, + } + + if verbose: + print( + f"[ok] field '{name}' written | " + f"shape={arr.shape} | p_shape={np.asarray(field.p).shape}" + ) + + return { + "comm": comm, + "decomp": decomp, + "fields": mugrid_fields, + "numpy": numpy_fields, + "report": report, + } + + +def print_field_summary(field_bundle: dict) -> None: + """ + Function that prints a summary of the muGrid fields stored in a field bundle. + + The printed summary includes the storage status and the shapes of the + non-ghost (`p`) and ghosted (`pg`) field views for each stored quantity. + + Parameters + ---------- + field_bundle : dict + Dictionary returned by `pack_adapted_grid_to_fields`. + + Returns + ------- + None + """ + print("\nmuGrid field summary") + print("-" * 72) + + for name, report in field_bundle["report"].items(): + print(f"{name}:") + print(f" status : {report['status']}") + print(f" shape : {report['shape']}") + print(f" p_shape : {report['p_shape']}") + print(f" pg_shape: {report['pg_shape']}") + print() + + +def pack_useful_fields_to_mugrid(result: dict, ghosts: int = 1, verbose: bool = False) -> dict: + """ + Function that stores only the most useful derived fields into muGrid containers. + + This is a reduced version of `pack_adapted_grid_to_fields` intended for + downstream workflows that only require a subset of the available fields. + + Parameters + ---------- + result : dict + Dictionary returned by `adapt_grid_to_circle`. + ghosts : int + Number of ghost points used in the muGrid decomposition. + verbose : bool + If True, print storage information during field creation. + + Returns + ------- + bundle : dict + Dictionary containing only the selected fields: + - "displacement" + - "inside" + - "interface" + - "dist" + - "k_node" + together with the corresponding communicator, decomposition, NumPy arrays, + and storage report. + """ + bundle = pack_adapted_grid_to_fields( + result, + ghosts=ghosts, + store_positions=False, + verbose=verbose, + ) + + useful_names = ["displacement", "inside", "interface", "dist", "k_node"] + + return { + "comm": bundle["comm"], + "decomp": bundle["decomp"], + "fields": {name: bundle["fields"][name] for name in useful_names}, + "numpy": {name: bundle["numpy"][name] for name in useful_names}, + "report": {name: bundle["report"][name] for name in useful_names}, + } + + +''' +//successful version + +""" +The code is genreated by Perplexity AI. + + +Reference: +Zecevic, M., Lebensohn, R. A., & Capolungo, L. (2026). +Achieving geometric accuracy in FFT-based micromechanical models +using conformal grid. Mechanics of Materials, 212, 105512. +[https://doi.org/10.1016/j.mechmat.2025.105512](https://doi.org/10.1016/j.mechmat.2025.105512) + + +""" + + +import numpy as np +from collections import deque +import muGrid + + + + +def make_grid_nodes(N: int = 64, L: float = 25.0) -> tuple[np.ndarray, np.ndarray, np.ndarray]: +    """ +    Construct a periodic-style Cartesian grid on [0, L) x [0, L). + + +    Notes +    ----- +    - The domain is interpreted as N cells per direction. +    - Only N x N grid points are stored. +    - The right and top boundary points are not stored explicitly. +    - Coordinates are stored with shape [xy, nx, ny]. + + +    Parameters +    ---------- +    N : int +        Number of cells per spatial direction. +    L : float +        Domain size. + + +    Returns +    ------- +    P : numpy.ndarray +        Grid-point coordinates with shape [xy, nx, ny]. +    x : numpy.ndarray +        x-coordinates of stored points, shape [nx]. +    y : numpy.ndarray +        y-coordinates of stored points, shape [ny]. +    """ +    x = np.linspace(0.0, L, N, endpoint=False) +    y = np.linspace(0.0, L, N, endpoint=False) +    X, Y = np.meshgrid(x, y, indexing="ij") +    P = np.stack([X, Y], axis=0) +    return P, x, y + + + + +def cell_labels( +    P: np.ndarray, +    center: tuple[float, float] = (0.0, 0.0), +    R: float = 20.0, +) -> np.ndarray: +    """ +    Label periodic cells as inside or outside a circular inclusion. + + +    Cell centers are built from the stored node at (i, j) and its periodic +    neighbors (i+1, j), (i, j+1), (i+1, j+1). + + +    Parameters +    ---------- +    P : numpy.ndarray +        Grid-point coordinates with shape [xy, nx, ny]. +    center : tuple of float +        Circle center. +    R : float +        Circle radius. + + +    Returns +    ------- +    inside : numpy.ndarray +        Integer array with shape [N, N]. +    """ +    cx, cy = center + + +    P_ip1 = np.roll(P, shift=-1, axis=1) +    P_jp1 = np.roll(P, shift=-1, axis=2) +    P_ip1_jp1 = np.roll(P_ip1, shift=-1, axis=2) + + +    P_ip1 = P_ip1.copy() +    P_jp1 = P_jp1.copy() +    P_ip1_jp1 = P_ip1_jp1.copy() + + +    P_ip1[0, -1, :] += R * 0.0 +    P_jp1[1, :, -1] += R * 0.0 + + +    P_ip1[0, -1, :] += (P[0, 1, 0] - P[0, 0, 0]) * P.shape[1] +    P_ip1_jp1[0, -1, :] += (P[0, 1, 0] - P[0, 0, 0]) * P.shape[1] +    P_jp1[1, :, -1] += (P[1, 0, 1] - P[1, 0, 0]) * P.shape[2] +    P_ip1_jp1[1, :, -1] += (P[1, 0, 1] - P[1, 0, 0]) * P.shape[2] + + +    Pc = 0.25 * (P + P_ip1 + P_jp1 + P_ip1_jp1) +    dx = Pc[0] - cx +    dy = Pc[1] - cy +    inside_mask = (dx * dx + dy * dy) <= R * R +    return inside_mask.astype(np.int32) + + + + +def interface_node_mask(cell_inside: np.ndarray) -> np.ndarray: +    """ +    Detect interface nodes on a periodic N x N stored grid. + + +    A stored point is marked if the surrounding periodic cells do not all have +    the same label. +    """ +    N = cell_inside.shape[0] +    mask = np.zeros((N, N), dtype=bool) + + +    for i in range(N): +        for j in range(N): +            vals = np.array([ +                cell_inside[(i - 1) % N, (j - 1) % N], +                cell_inside[i % N, (j - 1) % N], +                cell_inside[(i - 1) % N, j % N], +                cell_inside[i % N, j % N], +            ], dtype=int) +            if vals.min() != vals.max(): +                mask[i, j] = True + + +    return mask + + + + +def project_points_to_circle( +    Ppts: np.ndarray, +    center: tuple[float, float] = (0.0, 0.0), +    R: float = 20.0, +) -> np.ndarray: +    """ +    Project points onto the target circle. + + +    Accepts arrays with shape [xy, k] or [k, xy]. +    """ +    c = np.asarray(center, dtype=float) + + +    transposed = False +    if Ppts.ndim != 2: +        raise ValueError("Ppts must be a 2D array") +    if Ppts.shape[0] == 2: +        pts = Ppts.T +        transposed = True +    elif Ppts.shape[1] == 2: +        pts = Ppts +    else: +        raise ValueError("Ppts must have shape [xy, k] or [k, xy]") + + +    V = pts - c[None, :] +    r = np.linalg.norm(V, axis=1, keepdims=True) +    r = np.maximum(r, 1e-12) +    proj = c[None, :] + (R / r) * V +    return proj.T if transposed else proj + + + + +def manhattan_distance_to_interface(interface_mask: np.ndarray) -> np.ndarray: +    """ +    Compute periodic Manhattan distance to the nearest interface point. +    """ +    N = interface_mask.shape[0] +    dist = np.full((N, N), np.inf, dtype=float) +    queue = deque() + + +    for i in range(N): +        for j in range(N): +            if interface_mask[i, j]: +                dist[i, j] = 0.0 +                queue.append((i, j)) + + +    neighbors = [(1, 0), (-1, 0), (0, 1), (0, -1)] +    while queue: +        i, j = queue.popleft() +        for di, dj in neighbors: +            ii = (i + di) % N +            jj = (j + dj) % N +            if dist[ii, jj] > dist[i, j] + 1.0: +                dist[ii, jj] = dist[i, j] + 1.0 +                queue.append((ii, jj)) + + +    return dist + + + + +def stiffness_from_distance( +    dist: np.ndarray, +    k0: float = 1.0, +    b: float = 1.0, +    a: float = 1.0, +    kmin: float = 1e-3, +) -> np.ndarray: +    """ +    Compute nodal stiffness from distance to the interface. +    """ +    k = k0 / np.power(dist + a, b) +    return np.maximum(k, kmin) + + + + +def spring_relax_weighted( +    P: np.ndarray, +    fixed_mask: np.ndarray, +    k_node: np.ndarray, +    iters: int = 600, +    omega: float = 1.0, +) -> np.ndarray: +    """ +    Perform weighted spring relaxation on a periodic stored grid. + + +    Parameters +    ---------- +    P : numpy.ndarray +        Grid-point coordinates with shape [xy, nx, ny]. +    fixed_mask : numpy.ndarray +        Boolean mask with shape [nx, ny]. +    k_node : numpy.ndarray +        Stiffness field with shape [nx, ny]. +    iters : int +        Number of relaxation sweeps. +    omega : float +        Relaxation parameter. + + +    Returns +    ------- +    P_relaxed : numpy.ndarray +        Relaxed coordinates with shape [xy, nx, ny]. +    """ +    P_new = P.copy() +    N = P.shape[1] + + +    for _ in range(iters): +        for i in range(N): +            for j in range(N): +                if fixed_mask[i, j]: +                    continue + + +                neighbors = [ +                    ((i - 1) % N, j), +                    ((i + 1) % N, j), +                    (i, (j - 1) % N), +                    (i, (j + 1) % N), +                ] +                kij = k_node[i, j] +                w = np.empty(4, dtype=float) +                pts = np.empty((4, 2), dtype=float) + + +                xij = P_new[0, i, j] +                yij = P_new[1, i, j] + + +                for idx, (ii, jj) in enumerate(neighbors): +                    xnb = P_new[0, ii, jj] +                    ynb = P_new[1, ii, jj] + + +                    if ii == 0 and i == N - 1: +                        xnb += P[0, 1, 0] - P[0, 0, 0] +                        xnb += (N - 1) * (P[0, 1, 0] - P[0, 0, 0]) +                    elif ii == N - 1 and i == 0: +                        xnb -= P[0, 1, 0] - P[0, 0, 0] +                        xnb -= (N - 1) * (P[0, 1, 0] - P[0, 0, 0]) + + +                    if jj == 0 and j == N - 1: +                        ynb += P[1, 0, 1] - P[1, 0, 0] +                        ynb += (N - 1) * (P[1, 0, 1] - P[1, 0, 0]) +                    elif jj == N - 1 and j == 0: +                        ynb -= P[1, 0, 1] - P[1, 0, 0] +                        ynb -= (N - 1) * (P[1, 0, 1] - P[1, 0, 0]) + + +                    w[idx] = 0.5 * (kij + k_node[ii, jj]) +                    pts[idx] = [xnb, ynb] + + +                target = (w[:, None] * pts).sum(axis=0) / (w.sum() + 1e-15) +                P_new[:, i, j] = (1.0 - omega) * np.array([xij, yij]) + omega * target + + +    return P_new + + + + +def adapt_grid_to_circle( +    N: int = 64, +    L: float = 40.0, +    center: tuple[float, float] = (20.0, 20.0), +    R: float = 10.0, +    iters: int = 600, +    omega: float = 0.8, +    b: float = 1.0, +    k0: float = 0.25, +    kmin: float = 0.005, +) -> dict: +    """ +    Adapt a periodic-style grid to a circular inclusion. + + +    Notes +    ----- +    - There are N cells and only N x N stored points. +    - Right and top boundaries are not stored explicitly. +    - All coordinate arrays use shape [xy, nx, ny]. +    - Interface detection and distance are treated periodically. +    """ +    P0, _, _ = make_grid_nodes(N=N, L=L) +    inside = cell_labels(P0, center=center, R=R) +    interface = interface_node_mask(inside) +    fixed = interface.copy() + + +    P1 = P0.copy() +    idx = np.argwhere(interface) +    if idx.size: +        pts = P1[:, idx[:, 0], idx[:, 1]] +        P1[:, idx[:, 0], idx[:, 1]] = project_points_to_circle( +            pts, +            center=center, +            R=R, +        ) + + +    dist = manhattan_distance_to_interface(interface) +    k_node = stiffness_from_distance(dist, k0=k0, b=b, a=1.0, kmin=kmin) +    P2 = spring_relax_weighted(P1, fixed_mask=fixed, k_node=k_node, iters=iters, omega=omega) + + +    return { +        "P0": P0, +        "P": P2, +        "inside": inside, +        "interface": interface, +        "fixed": fixed, +        "dist": dist, +        "k_node": k_node, +        "params": { +            "N": N, +            "L": L, +            "center": center, +            "R": R, +            "iters": iters, +            "omega": omega, +            "b": b, +            "k0": k0, +            "kmin": kmin, +        }, +    } + + + +if __name__ == "__main__": +    coor, x, y = make_grid_nodes(N=4, L=1.0) +    print(coor[:, 0, 0]) +    print(coor[:, 1, 0]) + + +    nx = 8 +    ny = 8 +    g = 1 +    print("before communicator") +    comm = muGrid.Communicator() +    print("after communicator") +    print("before decomp") + + +    decomp = muGrid.CartesianDecomposition( +        communicator=comm, +        nb_domain_grid_pts=(nx, ny), +        nb_subdivisions=(1, 1), +        nb_ghosts_left=(g, g), +        nb_ghosts_right=(g, g), +    ) +    print("after decomp") + + +    field = decomp.real_field("displacement", components=(3,)) +    print(comm) +    print(field) +    print(type(field)) +    print([name for name in dir(field) if not name.startswith("_")]) + + +''' From 44c735938ff3d106ae302936b0974c270b033bec Mon Sep 17 00:00:00 2001 From: MartinLadecky Date: Mon, 6 Jul 2026 14:26:01 +0200 Subject: [PATCH 2/3] adding new experiment with multiple phases --- ...ty_transformed_grid_Jia_multiple_phases.py | 217 ++++++++++++++++++ 1 file changed, 217 insertions(+) create mode 100644 experiments/grid_adaptation/exp_2D_homogenization_conductivity_transformed_grid_Jia_multiple_phases.py 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}') From e14649886555994151689b55a7a7839776231af9 Mon Sep 17 00:00:00 2001 From: MartinLadecky Date: Wed, 16 Sep 2026 18:23:28 +0200 Subject: [PATCH 3/3] add otsu fo gia --- muFFTTO/otsu.py | 421 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 421 insertions(+) create mode 100644 muFFTTO/otsu.py 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