Make a script to print out the information of all bubbles

%% Cell type:code id:600f8b20-7dcf-40d3-82f2-011349965da3 tags:

``` python

# Some parameters for the behavior of the extraction

# Specifies the spacing of the sampling grid.

lattice_parameter = 3.177

# The radius to use for a He splat

He_radius = 0.01

# Minimum number of He atoms that have to escape for a bubble to be

# considered burst.

burst_size = 1

# Minimum fraction of He atoms in a bubble that have to escape for

# the buble to be considered burst.

burst_fraction = 0.75

# Number of time steps to pass to check for escaped atoms.

time_to_escape = 10

# The maximum distance (squared) between bubble centers in two

# different timesteps for them to be considered the same bubble.

same_distance = 16

# Directory containing VTK files holding atoms.

#data_path = '/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm/vtk-files'

data_path = '/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files'

```

%% Cell type:code id:d2407a13-cf90-4203-92cb-c716928c96a4 tags:

``` python

import sys

sys.path.append('/Users/4d5/builds/VTK/Release/lib/python3.11/site-packages')

import vtk

import numpy as np

from vtk.util.numpy_support import vtk_to_numpy

```

%% Cell type:code id:67821616-203b-4523-8459-35c4e2700a2d tags:

``` python

class Bubble:

    '''Holds data for a bubble of atoms.

Bubble is constructed with a vtkDataSet that contains the atoms

of the bubble. Alternately, it can be constructed with a set of

classified atoms and the identifier of the atom in the `RegionId`

field (as comes from find_bubbles).'''

    def __init__(self, data, bubble_id = 0):

        if (bubble_id < 1):

            self.data = data

        else:

            threshold = vtk.vtkThreshold()

            threshold.SetInputData(data)

            threshold.SetInputArrayToProcess(

                0, 0, 0,

                'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'RegionId')

            threshold.SetLowerThreshold(bubble_id - 0.1)

            threshold.SetUpperThreshold(bubble_id + 0.1)

            threshold.Update()

            self.data = threshold.GetOutput()

        self._center = None

    def get_number_of_atoms(self):

        return self.data.GetNumberOfPoints()

    def get_center(self):

        if not self._center:

            compute_center = vtk.vtkCenterOfMass()

            compute_center.SetInputData(self.data)

            compute_center.Update()

            self._center = compute_center.GetCenter()

        return self._center

```

%% Cell type:code id:364b68d8-1c5f-4e00-af73-15242a059e10 tags:

``` python

class HeAtoms:

    '''Holds information about one timestep of He atoms.

This class is initialized by providing a VTK dataset representing

a LAMMPS atom file containing an `id` field giving unique identifiers

to atoms and a `type` field identifying the type of each atom

(`type` == 2 is He atoms). The data needs no cell data. The constructor

must also be given a timestep.

Use the `ReadVTKFile()` method to load the data from a VTK data set.

The `bubbles` field contains a list of `Bubble` objects holding the

extracted bubbles. The `timestep` field contains the time of the

atoms. The `filename` field provides the file the atoms came from,

when available.'''

    def ReadVTKFile(filename):

        import re

        timestep_match = re.search(r'([0-9]+)\.vtk', filename)

        timestep = int(timestep_match.group(1))

        reader = vtk.vtkDataSetReader()

        reader.SetFileName(filename)

        reader.Update()

        return HeAtoms(reader.GetOutput(), timestep, filename)

    def __init__(self, data, timestep, filename=None):

        self.timestep = timestep

        self.filename = filename

        # Many filters expect the data to have cells. Make one cell per point.

        tocloud = vtk.vtkConvertToPointCloud()

        tocloud.SetInputData(data)

        tocloud.SetCellGenerationMode(tocloud.VERTEX_CELLS)

        # Extract the He atoms.

        extract_he = vtk.vtkThreshold()

        extract_he.SetInputConnection(tocloud.GetOutputPort())

        extract_he.SetInputArrayToProcess(

            0, 0, 0, 'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'type')

        extract_he.SetLowerThreshold(1.9)

        extract_he.SetUpperThreshold(2.1)

        # Create a density field of atoms by splatting the atoms into a grid.

        # Here we use `vtkGaussianSplatter`. If this is taking too long, it

        # is possible to speed things up using a `vtkFastSplatter`.

        bounds = data.GetBounds()

        sampling_resolution = (

            int((bounds[1] - bounds[0]) / (lattice_parameter / 2)),

            int((bounds[3] - bounds[2]) / (lattice_parameter / 2)),

            int((bounds[5] - bounds[4]) / (lattice_parameter / 2)))

        splatter = vtk.vtkGaussianSplatter()

        splatter.SetInputConnection(extract_he.GetOutputPort())

        splatter.SetSampleDimensions(sampling_resolution)

        splatter.SetRadius(He_radius)

        splatter.SetNormalWarping(0)

        splatter.SetAccumulationModeToSum()

        # Use an image connectivity filter to find regions with nearby atoms.

        # The filter will assign each connected region a unique ID and write

        # that to point data.

        connectivity = vtk.vtkImageConnectivityFilter()

        connectivity.SetInputConnection(splatter.GetOutputPort())

        connectivity.SetLabelScalarTypeToInt()

        #connectivity.SetLabelModeToSizeRank()

        # Lowering threshold to include regions to prevent missing atoms at edge

        # of bubble.

        connectivity.SetScalarRange(0.1, 100000000000)

        # Frustratingly, there appears to be a bug in `vtkImageConnectivityFilter`

        # that does not properly pass the origin and spacing of the image, which

        # means that it will no longer align with the atom data. Restore it to

        # the proper values.

        splatter.Update()

        splatter_output = splatter.GetOutput()

        reset_position = vtk.vtkImageChangeInformation()

        reset_position.SetInputConnection(connectivity.GetOutputPort())

        reset_position.SetOutputOrigin(splatter_output.GetOrigin())

        reset_position.SetOutputSpacing(splatter_output.GetSpacing())

        # Sample the connectivity field back onto the atoms that generated them.

        # This will add a `RegionId` field to the atoms that identifes which

        # cluster each one belongs to.

        probe = vtk.vtkProbeFilter()

        probe.SetInputConnection(0, extract_he.GetOutputPort())

        probe.SetInputConnection(1, reset_position.GetOutputPort())

        probe.CategoricalDataOn()

        probe.PassPointArraysOn()

        probe.Update()

        self.atoms = probe.GetOutput()

        self.bubbles = []

        for bubble_index in range(1, int(self.atoms.GetPointData().

                                         GetArray('RegionId').GetRange()[1])):

            bubble = Bubble(self.atoms, bubble_index)

            if bubble.get_number_of_atoms() >= burst_size:

                self.bubbles.append(bubble)

```

%% Cell type:code id:02aeb93b-6cb3-448c-9911-c2a0df9834e5 tags:

``` python

timestep = '084001'

sorted_atoms = HeAtoms.ReadVTKFile(data_path + '/110-25x25nm-detailed-' + timestep + '.vtk')

```

%% Cell type:code id:66149a82-9735-43a6-8a0e-375888a60b51 tags:

``` python

f = open('bubbles-' + timestep + '.csv', 'wt', encoding='utf-8')

f.write('center_x,center_y,center_z,size\n')

for bubble in sorted_atoms.bubbles:

    center = bubble.get_center()

    f.write('{},{},{},{}\n'.format(

        center[0], center[1], center[2], bubble.get_number_of_atoms()))

f.close()

```

%% Cell type:code id:3da1430f-dde9-4cf8-8b4b-cebdce2c9b8b tags:

``` python

```