diff --git a/cipher_parse/cipher_input.py b/cipher_parse/cipher_input.py index a7d5a4d..9639dbc 100644 --- a/cipher_parse/cipher_input.py +++ b/cipher_parse/cipher_input.py @@ -333,6 +333,7 @@ def from_voronoi( num_phases=None, random_seed=None, is_periodic=False, + combine_phases=None, ): geometry = CIPHERGeometry.from_voronoi( num_phases=num_phases, @@ -343,6 +344,7 @@ def from_voronoi( size=size, random_seed=random_seed, is_periodic=is_periodic, + combine_phases=combine_phases, ) inp = cls( @@ -366,6 +368,7 @@ def from_seed_voronoi( solution_parameters, random_seed=None, is_periodic=False, + combine_phases=None, ): return cls.from_voronoi( seeds=seeds, @@ -378,6 +381,7 @@ def from_seed_voronoi( solution_parameters=solution_parameters, random_seed=random_seed, is_periodic=is_periodic, + combine_phases=combine_phases, ) @classmethod @@ -393,6 +397,7 @@ def from_random_voronoi( solution_parameters, random_seed=None, is_periodic=False, + combine_phases=None, ): return cls.from_voronoi( num_phases=num_phases, @@ -405,6 +410,7 @@ def from_random_voronoi( solution_parameters=solution_parameters, random_seed=random_seed, is_periodic=is_periodic, + combine_phases=combine_phases, ) @classmethod @@ -418,6 +424,7 @@ def from_voxel_phase_map( outputs, solution_parameters, random_seed=None, + combine_phases=None, ): geometry = CIPHERGeometry( voxel_phase=voxel_phase, @@ -425,6 +432,7 @@ def from_voxel_phase_map( interfaces=interfaces, size=size, random_seed=random_seed, + combine_phases=combine_phases, ) inp = cls( geometry=geometry, @@ -445,6 +453,7 @@ def from_dream3D( solution_parameters, container_labels=None, phase_type_map=None, + combine_phases=None, ): default_container_labels = { "SyntheticVolumeDataContainer": "SyntheticVolumeDataContainer", @@ -559,6 +568,7 @@ def from_dream3D( components=components, outputs=outputs, solution_parameters=solution_parameters, + combine_phases=combine_phases, ) @property diff --git a/cipher_parse/geometry.py b/cipher_parse/geometry.py index 79df855..e260547 100644 --- a/cipher_parse/geometry.py +++ b/cipher_parse/geometry.py @@ -1,3 +1,5 @@ +import random + from damask import Orientation import pyvista as pv import numpy as np @@ -39,6 +41,7 @@ def __init__( time=None, increment=None, incremental_data_idx=None, + combine_phases=None, ): """ Parameters @@ -107,6 +110,22 @@ def __init__( self._ensure_phase_assignment(random_seed) + if combine_phases: + # combine multiple phases into the same phase to reduce the computational cost + self.voxel_phase = self.combine_phases_per_phase_type( + voxel_map, + materials, + combine_phases, + ) + self.voxel_map = VoxelMap( + region_ID=voxel_phase, + size=size, + is_periodic=is_periodic, + quiet=quiet, + ) + all_phases = self.present_phases + self._num_phases = all_phases.size + self._phase_material = self._get_phase_material() self._validate_interfaces() self._check_interface_phase_pairs() @@ -131,6 +150,93 @@ def __init__( self._misorientation_matrix = None self._misorientation_matrix_is_degrees = None + @staticmethod + def combine_phases_per_phase_type(voxel_map, materials, combine_phases): + + print(f"combining phases according to {combine_phases}") + + voxel_phase_new = voxel_map.region_ID + neighbours = voxel_map.neighbour_list + + for mat in materials: + for pt_i in mat.phase_types: + print(f"{pt_i.name=}") + max_phases_i = combine_phases.get(pt_i.name) + if max_phases_i: + + n_phases = pt_i.phases.size + group_size = int(np.ceil(n_phases / max_phases_i)) + + print(f" n_phases: {n_phases}") + print(f" max_phases_i: {max_phases_i}") + print(f" group_size: {group_size}") + + voxel_phase_new = np.copy(voxel_phase_new) + + reassignment = {} + reassignment_inv = {} + possible_IDs_idx = set(range(max_phases_i, n_phases)) + possible_IDs = set(pt_i.phases[list(possible_IDs_idx)]) + count = 0 + kept_IDs = [] + for root_ID_idx in range(int(np.ceil(n_phases / group_size))): + + root_ID = pt_i.phases[root_ID_idx] + kept_IDs.append(root_ID) + + neighbours_i = set(neighbours[:, neighbours[0] == root_ID][1]) + possible_IDs -= {root_ID} + possible_IDs_i = possible_IDs - neighbours_i + shared_IDs = [root_ID] + count += 1 + + for _ in range(1, group_size): + if count >= n_phases: + break + try: + sampled_ID = random.sample(possible_IDs_i, 1)[0] + except ValueError: + print( + f"No non-neighbouring samples left for root_ID: {root_ID}." + ) + # allow touching phase IDs within this group: + possible_IDs_i = possible_IDs - set(shared_IDs) + + if possible_IDs_i: + sampled_ID = random.sample(possible_IDs_i, 1)[0] + else: + raise ValueError( + f"Cannot find non-neighbouring root IDs" + ) from None + + neighbours_sampled = set( + neighbours[:, neighbours[0] == sampled_ID][1] + ) + possible_IDs_i -= neighbours_sampled + possible_IDs_i -= {sampled_ID} + possible_IDs -= {sampled_ID} + shared_IDs.append(sampled_ID) + count += 1 + + for i in shared_IDs: + reassignment[i] = root_ID + voxel_phase_new[voxel_phase_new == i] = root_ID + + reassignment_inv[root_ID] = shared_IDs + + pt_i.phases = kept_IDs # modify phase type phases + + # reindex phases across all materials to maintain consecutive phase IDs: + voxel_phase_new_flat = voxel_phase_new.reshape(-1) + uniq, inv = np.unique(voxel_phase_new_flat, return_inverse=True) + reindex = dict(zip(uniq, range(len(uniq)))) + voxel_phase_new = inv.reshape(voxel_phase_new.shape) + for mat_j in materials: + for pt_j in mat_j.phase_types: + pt_j.phases = np.array([reindex[i] for i in pt_j.phases]) + + return voxel_phase_new + def __eq__(self, other): # Note we don't check seeds (not stored in YAML file) if not isinstance(other, self.__class__): @@ -456,7 +562,7 @@ def _get_phase_material(self): "Not all phases are accounted for in the phase type definitions." ) # TODO: test raise - # check all phase indices form a consequtive range: + # check all phase indices form a consecutive range: num_phases_range = set(np.arange(self.num_phases)) known_phases = set(np.hstack(all_phase_idx)) miss_phase_idx = num_phases_range - known_phases @@ -740,6 +846,7 @@ def from_voronoi( num_phases=None, random_seed=None, is_periodic=False, + combine_phases=None, ): if sum(i is not None for i in (seeds, num_phases)) != 1: raise ValueError(f"Specify exactly one of `seeds` and `num_phases`") @@ -769,6 +876,7 @@ def from_voronoi( size=size, seeds=seeds, random_seed=random_seed, + combine_phases=combine_phases, ) @classmethod