Multistep continuous simulation with regridding

[mehady007]

Reference: Redrid
Regrids a DAMASK simulation by mapping the final state of an existing result onto a new grid and generating compatible geometry and restart files for seamless continuation.

import re
import numpy as np
import h5py
import damask
from pathlib import Path
rng = np.random.default_rng()

class RegridPreproces:
    """
    Pre-processing utility for regridding DAMASK simulations.

    This class enables continuation of a DAMASK simulation on a new grid
    resolution by using:
      - a completed result file (.hdf5),
      - the corresponding restart file,
      - and the original geometry (.vti).

    The workflow:
      1. Extracts deformation from the final increment of a result file.
      2. Determines a new grid resolution (automatic, scaled, or manual).
      3. Computes a grid-to-grid mapping using DAMASK regridding filters.
      4. Generates a new geometry file consistent with the deformed state.
      5. Constructs a new restart file with all state variables mapped
         onto the new grid.

    This allows simulations to continue seamlessly after mesh refinement,
    resolution changes, or large deformations without restarting from
    the initial configuration.
    """
    def __init__(self,
                 data_folder=None,
                 ):

        self.data_folder = Path(data_folder) if data_folder else Path.cwd()


        files = sorted(f for f in self.data_folder.glob("*.hdf5")
                       if not f.name.endswith("_restart.hdf5"))

        if not files:
            raise FileNotFoundError("No .hdf5 result files found.")

        print("\n[📂] Available result files:")
        for i, f in enumerate(files, 1):
            print(f"  {i}: {f.name}")
        try:
            idx = int(input(f"Select result file [1-{len(files)}]: "))
            assert 1 <= idx <= len(files)
        except Exception:
            raise ValueError("Invalid selection.")
        self.result_file = files[idx-1]
        print(f"[✔] Selected: {self.result_file.name}")
        

    def ensure_unique_filename(self,grid_in, grid_out):
        if grid_in != grid_out:
            return grid_in, grid_out        
        pass_num = self.pass_num
         # strip any existing _vN before re-appending
        grid_old_base = re.sub(r"_v\d+$", "", grid_in)
        grid_new_base = re.sub(r"_v\d+$", "", grid_out)        
        if pass_num <=1:

            grid_old = f"{grid_old_base}"
            grid_new = f"{grid_new_base}_v{pass_num + 1}"           
        else:   
            grid_old = f"{grid_old_base}_v{pass_num}"
            grid_new = f"{grid_new_base}_v{pass_num + 1}"

        return grid_old, grid_new
    
    def generate_mapping_from_result(self, scale_vec=None, cells_manual=None):
        stem = self.result_file.stem  # e.g. something_step1 or pass1
        m = re.match(r"(\d+)(.*?)(\d+x\d+x\d+)_seed(\d+)(.*?)(\d+)$", stem)
        if not m:
            raise ValueError("Result filename must end with a number (e.g., step1.hdf5).")
        self.N = int(m.group(1))        
        self.seed = int(m.group(4))
        self.pass_num = int(m.group(6))
        self.output_number = self.pass_num  + 1
        self.fname_in  = self.data_folder / f"{m.group(0)}_restart.hdf5"

        print(f"[📂] Loading result file: {self.result_file.name}")
        r = damask.Result(str(self.result_file))
        initial_cells = np.array(r.cells, dtype=int)
        old_vti_file = f"{self.N}{m.group(2)}{m.group(3)}_seed{self.seed}.vti"

        F_last = r.view(increments=-1).get('F')  # (Nx,Ny,Nz,3,3)
        F_avg  = np.average(F_last, axis=0)

        if cells_manual is not None:
            cells_new = np.array(cells_manual, dtype=int)
            print(f"[🔢] Using MANUAL cells: {cells_new}")
        elif scale_vec is not None:
            if len(scale_vec) != 3:
                raise ValueError("scale_vec must be length 3 [x,y,z].")
            cells_new = np.ceil(initial_cells * np.array(scale_vec)).astype(int)
            print(f"[🔢] Using user scaling {scale_vec} => new cells: {cells_new}")
        else:
            diag = np.diag(F_avg)
            cells_new = np.ceil(initial_cells * diag).astype(int)
            print(f" Using ⟨F⟩ diag={diag} => new cells: {cells_new}")

        if cells_new[2] % 2 != 0:
            cells_new[2] += 1
        self.cells_new = cells_new
        self.grid_old = f"{self.N}{m.group(2)}{m.group(3)}_seed{self.seed}"
        self.grid_new = f"{self.N}{m.group(2)}{cells_new[0]}x{cells_new[1]}x{cells_new[2]}_seed{self.seed}"
        self.fname_out = self.data_folder / f"{self.grid_new}{m.group(5)}{self.output_number}_restart.hdf5"
        self.grid_old, self.grid_new = self.ensure_unique_filename(self.grid_old, self.grid_new)
       
        print(f" Geometry file flow: {self.grid_old}.vti → {self.grid_new}.vti")
        print(f" Restart file flow: \n -> {getattr(self.fname_in, 'name')}\n -> {getattr(self.fname_out, 'name')}")            
        # compute mapping
        mapping = damask.grid_filters.regrid(
            r.size,
            F_last.reshape(tuple(initial_cells) + (3,3)),
            self.cells_new
        )
        self.mapping_flat = mapping.flatten(order='F')
        

        # save new VTI
        g_old = damask.GeomGrid.load(str(self.data_folder / f"{self.grid_old}.vti"))
        g_old.size = (F_avg*np.eye(3)) @ g_old.size
        g_old.assemble(mapping).save(str(self.data_folder / f"{self.grid_new}.vti"))
        print(f"[💾] Regridded geometry saved: {self.grid_new}.vti")
        return mapping    
    def regrid_restart(self):
        with (h5py.File(self.fname_in) as f_in, h5py.File(self.fname_out,'w') as f_out):
            f_in.copy('homogenization',f_out)
    
            f_out.create_group('phase')
            for label in f_in['phase']:
                f_out['phase'].create_group(label)
                F_p = damask.tensor.transpose(f_in['phase'][label]['F_p'][()])
                F = damask.tensor.transpose(f_in['phase'][label]['F'][()])
                f_out['phase'][label].create_dataset('F',data=np.broadcast_to(np.eye(3),(len(self.mapping_flat),3,3,)))
                f_out['phase'][label].create_dataset('F_p',data=damask.tensor.transpose(np.matmul(F_p,np.linalg.inv(F)))[self.mapping_flat])
                for d in f_in['phase'][label]:
                    if d in f_out[f'phase/{label}']: continue
                    f_out['phase'][label].create_dataset(d,data=f_in['phase'][label][d][()][self.mapping_flat])
    
            f_out.create_group('solver')
            for d in ['F','F_lastInc']:
                f_out['solver'].create_dataset(d,data=np.broadcast_to(np.eye(3),np.append(self.cells_new.prod(),(3,3))))
            for d in ['F_aim', 'F_aim_lastInc']:
                f_out['solver'].create_dataset(d,data=np.eye(3))
            f_out['solver'].create_dataset('F_aimDot',data=np.zeros((3,3)))
            for d in f_in['solver']:
                if d not in f_out['solver']: f_in['solver'].copy(d,f_out['solver'])    

How to use:

dat_folder= Path('your working path')
obj = RegridPreproces(dat_folder)
mapping = obj.generate_mapping_from_result()
obj.regrid_restart()
np.savetxt(dat_folder/'map.txt',mapping.flatten('F'),fmt='%d')