3D foam example

[Kresa96]

Hello,

I like to build a structure very similar to the foam example shown on the documentation page (see image). How to modify the corresponding code to built this kind of structures in 3D? It would be also very nice to implement a feature, where one can simply put spherical particles into a matrix of another phase and translate it directly into a voxel grid/ numpy array without the voronoi step.

[kip-hart]

Hi @Kresa96,

You could avoid the voronoi step by creating a list of void spheres, then test whether each voxel is inside a sphere or not.

For example, the following code would create the voids:

import microstructpy as msp
import numpy as np
import scipy.stats

# Create Voids
material_1 = {
    'name': 'Voids',
    'material_type': 'void',
    'fraction': 1,
    'shape': 'sphere',
    'size': scipy.stats.uniform(loc=0, scale=1.5)
}
# Create Domain
domain = msp.geometry.Cube(side_length=15, corner=(0, 0))

# Create List of Un-Positioned Seeds
porosity = 0.12;
seed_volume = porosity * domain.volume
rng_seeds = {'size': 1}
seeds = msp.seeding.SeedList.from_info([material_1],
                                       seed_volume,
                                       rng_seeds)

# Position Seeds in Domain
seeds.position(domain)

Then to determine which voxels are within a void:

# 1. Create node and element grids
mesh_size = 1
mins, maxs = zip(*domain.limits)
sides = [np.arange(lb, ub, mesh_size) for lb, ub in zip(mins, maxs)]
mgrid = np.meshgrid(*sides)
nodes = np.array([g.flatten() for g in mgrid]).T
node_nums = np.arange(mgrid[0].size).reshape(mgrid[0].shape)

m, n, p = node_nums.shape
kp1 = node_nums[:(m-1), :(n-1), :(p-1)].flatten()
kp2 = node_nums[1:m, :(n-1), :(p-1)].flatten()
kp3 = node_nums[1:m, 1:n, :(p-1)].flatten()
kp4 = node_nums[:(m-1), 1:n, :(p-1)].flatten()
kp5 = node_nums[:(m-1), :(n-1), 1:p].flatten()
kp6 = node_nums[1:m, :(n-1), 1:p].flatten()
kp7 = node_nums[1:m, 1:n, 1:p].flatten()
kp8 = node_nums[:(m-1), 1:n, 1:p].flatten()
elems = np.array([kp1, kp2, kp3, kp4, kp5, kp6, kp7, kp8]).T

# 2. Compute element centers
cens = nodes[elems[:, 0]] + 0.5 * mesh_size

# 3. For each seed:
i_remain = np.arange(cens.shape[0])
elem_regs = np.full(cens.shape[0], -1)
elem_atts = np.full(cens.shape[0], -1)
for s_num, seed in enumerate(seeds):
    # A. Create a bounding box
    [s_mins, s_maxs] = zip(*seed.geometry.limits)

    # B. Isolate element centers with box
    s_i_remain = np.copy(i_remain)
    for i, lb in enumerate(s_mins):
        ub = s_maxs[i]
        x = cens[s_i_remain, i]
        in_range = (x >= lb) & (x <= ub)
        s_i_remain = s_i_remain[in_range]
    
    # C. Assign remaining centers to seed
    elem_regs[s_i_remain] = s_num
    elem_atts[s_i_remain] = s_num
    i_remain = np.setdiff1d(i_remain, s_i_remain)

# 4. Remove voids
# Remove element
rm_mask = elem_atts == -1
elems = elems[rm_mask]
elem_atts = elem_atts[rm_mask]

# Re-number nodes
nodes_mask = np.isin(np.arange(nodes.shape[0]), elems)
n_remain = np.sum(nodes_mask)
node_n_conv = np.arange(nodes.shape[0])
node_n_conv[nodes_mask] = np.arange(n_remain)

nodes = nodes[nodes_mask]
elems = node_n_conv[elems]

You can then use these nodes and elements to create a RasterMesh object and take advantage of its methods.

rmesh = msp.meshing.RasterMesh(nodes, elems)

One word of caution with foam meshes is that it's very easy to create disconnected regions of material. In the example you show from the documentation, there are extra steps involved to help prevent that. You may need to play around with the scripts to make it consistently generate voxel meshes without disconnected regions.

Let me know if I can support further.