Source code for sfepy.homogenization.recovery

from __future__ import print_function
from __future__ import absolute_import
import os

import numpy as nm

from sfepy.base.base import get_default, Struct, output
from sfepy.base.ioutils import get_print_info
from sfepy.base.timing import Timer
from sfepy.discrete.fem import extend_cell_data, Mesh
from sfepy.homogenization.utils import coor_to_sym
from sfepy.base.conf import get_standard_keywords
from sfepy.discrete import Problem, Region
from sfepy.base.conf import ProblemConf
from sfepy.homogenization.coefficients import Coefficients
from sfepy.homogenization.micmac import get_correctors_from_file_hdf5
import os.path as op
import six
from six.moves import range
import atexit

shared = Struct()

#
# TODO : interpolate fvars to macro times. ?mid-points?
#
# TODO : clean-up!
#


[docs] def get_output_suffix(iel, ts, naming_scheme, format, output_format): if output_format != 'h5': if naming_scheme == 'step_iel': suffix = '.'.join((ts.suffix % ts.step, format % iel)) else: suffix = '.'.join((format % iel, ts.suffix % ts.step)) else: suffix = format % iel return suffix
[docs] def convolve_field_scalar(fvars, pvars, iel, ts): r""" .. math:: \int_0^t f(t-s) p(s) ds Notes ----- - t is given by step - f: fvars scalar field variables, defined in a micro domain, have shape [step][fmf dims] - p: pvars scalar point variables, a scalar in a point of macro-domain, FMField style have shape [n_step][var dims] """ step0 = max(0, ts.step - fvars.steps[-1]) val = nm.zeros_like(fvars[0]) for ik in range(step0, ts.step + 1): vf = fvars[ts.step-ik] vp = pvars[ik][iel, 0, 0, 0] val += vf * vp * ts.dt return val
[docs] def convolve_field_sym_tensor(fvars, pvars, var_name, dim, iel, ts): r""" .. math:: \int_0^t f^{ij}(t-s) p_{ij}(s) ds Notes ----- - t is given by step - f: fvars field variables, defined in a micro domain, have shape [step][fmf dims] - p: pvars sym. tensor point variables, a scalar in a point of macro-domain, FMField style, have shape [dim, dim][var_name][n_step][var dims] """ step0 = max(0, ts.step - fvars[0, 0][var_name].steps[-1]) val = nm.zeros_like(fvars[0, 0][var_name][0]) for ik in range(step0, ts.step + 1): for ir in range(dim): for ic in range(dim): ii = coor_to_sym(ir, ic, dim) vf = fvars[ir, ic][var_name][ts.step-ik] vp = pvars[ik][iel, 0, ii, 0] val += vf * vp * ts.dt return val
[docs] def add_strain_rs(corrs_rs, strain, vu, dim, iel, out=None): if out is None: out = nm.zeros_like(corrs_rs[0, 0][vu][0]) for ir in range(dim): for ic in range(dim): ii = coor_to_sym(ir, ic, dim) out += corrs_rs[ir, ic][vu].data * strain[iel, 0, ii, 0] return out
[docs] def combine_scalar_grad(corrs, grad, vn, ii, shift_coors=None): r""" .. math:: \eta_k \partial_k^x p or .. math:: (y_k + \eta_k) \partial_k^x p """ dim = grad.shape[2] if shift_coors is None: out = corrs[0][vn].data * grad[ii, 0, 0, 0] for ir in range(1, dim): out += corrs[ir][vn].data * grad[ii, 0, ir, 0] else: out = (shift_coors[:, 0] + corrs[0][vn].data) * grad[ii, 0, 0, 0] for ir in range(1, dim): out += (shift_coors[:, ir] + corrs[ir][vn].data) \ * grad[ii, 0, ir, 0] return out
[docs] def compute_u_corr_steady(corrs_rs, strain, vu, dim, iel): r""" .. math:: \sum_{ij} \bm{\omega}^{ij}\, e_{ij}(\bm{u}) Notes ----- - iel = element number """ u_corr = add_strain_rs(corrs_rs, strain, vu, dim, iel) return u_corr
[docs] def compute_u_corr_time(corrs_rs, dstrains, corrs_pressure, pressures, vu, dim, iel, ts): r""" .. math:: \sum_{ij} \left[ \int_0^t \bm{\omega}^{ij}(t-s) {\mathrm{d} \over \mathrm{d} s} e_{ij}(\bm{u}(s))\,ds\right] + \int_0^t \widetilde{\bm{\omega}}^P(t-s)\,p(s)\,ds """ u_corr = convolve_field_scalar(corrs_pressure[vu], pressures, iel, ts) u_corr += convolve_field_sym_tensor(corrs_rs, dstrains, vu, dim, iel, ts) return u_corr
[docs] def compute_p_corr_steady(corrs_pressure, pressure, vp, iel): r""" .. math:: \widetilde\pi^P\,p """ p_corr = corrs_pressure[vp].data * pressure[iel, 0, 0, 0] return p_corr
[docs] def compute_p_corr_time(corrs_rs, dstrains, corrs_pressure, pressures, vdp, dim, iel, ts): r""" .. math:: \sum_{ij} \int_0^t {\mathrm{d} \over \mathrm{d} t} \widetilde\pi^{ij}(t-s)\, {\mathrm{d} \over \mathrm{d} s} e_{ij}(\bm{u}(s))\,ds + \int_0^t {\mathrm{d} \over \mathrm{d} t}\widetilde\pi^P(t-s)\,p(s)\,ds """ p_corr = convolve_field_scalar(corrs_pressure[vdp], pressures, iel, ts) p_corr += convolve_field_sym_tensor(corrs_rs, dstrains, vdp, dim, iel, ts) return p_corr
[docs] def compute_u_from_macro(strain, coor, iel, centre=None): r""" Macro-induced displacements. .. math:: e_{ij}^x(\bm{u})\,(y_j - y_j^c) """ n_nod, dim = coor.shape if centre is None: centre = nm.zeros((dim,), dtype=nm.float64) n_nod, dim = coor.shape um = nm.zeros((n_nod * dim,), dtype=nm.float64) for ir in range(dim): for ic in range(dim): ii = coor_to_sym(ir, ic, dim) um[ir::dim] += strain[iel, 0, ii, 0] * (coor[:, ic] - centre[ic]) return um
[docs] def compute_p_from_macro(p_grad, coor, iel, centre=None, extdim=0): r""" Macro-induced pressure. .. math:: \partial_j^x p\,(y_j - y_j^c) """ n_nod, dim = coor.shape if centre is None: centre = nm.zeros((dim,), dtype=nm.float64) n_nod, dim = coor.shape pm = nm.zeros((n_nod,), dtype=nm.float64) for ic in range(dim + extdim): pm += p_grad[iel, 0, ic, 0] * (coor[:, ic] - centre[ic]) return pm
[docs] def compute_micro_u(corrs, strain, vu, dim, out=None): r""" Micro displacements. .. math:: \bm{u}^1 = \bm{\chi}^{ij}\, e_{ij}^x(\bm{u}^0) """ if out is None: out = nm.zeros_like(corrs[vu+'_00']) for ir in range(dim): for ic in range(dim): ii = coor_to_sym(ir, ic, dim) out += corrs[vu+'_%d%d' % (ir, ic)] * strain[ii] return out
[docs] def compute_stress_strain_u(pb, integral, region, material, vu, data): var = pb.create_variables([vu])[vu] var.set_data(data) stress = pb.evaluate('ev_cauchy_stress.%s.%s(%s, %s)' % (integral, region, material, vu), verbose=False, mode='el_avg', **{vu: var}) strain = pb.evaluate('ev_cauchy_strain.%s.%s(%s)' % (integral, region, vu), verbose=False, mode='el_avg', **{vu: var}) return extend_cell_data(stress, pb.domain, region), \ extend_cell_data(strain, pb.domain, region)
[docs] def add_stress_p(out, pb, integral, region, vp, data): var = pb.create_variables([vp])[vp] var.set_data(data) press0 = pb.evaluate('ev_integrate.%s.%s(%s)' % (integral, region, vp), verbose=False, mode='el_avg', **{vp: var}) press = extend_cell_data(press0, pb.domain, region) dim = pb.domain.mesh.dim nn = out.shape[0] for ii in range(nn): for j in range(dim): out[ii, 0, j, 0] += press[ii, 0, 0, 0]
[docs] def compute_mac_stress_part(pb, integral, region, material, vu, mac_strain): avgmat = pb.evaluate('ev_integrate_mat.%s.%s(%s, %s)' % (integral, region, material, vu), verbose=False, mode='el_avg') return extend_cell_data(nm.dot(avgmat, mac_strain), pb.domain, region)
[docs] def recover_bones(problem, micro_problem, region, eps0, ts, strain, dstrains, p_grad, pressures, corrs_permeability, corrs_rs, corrs_time_rs, corrs_pressure, corrs_time_pressure, var_names, naming_scheme='step_iel'): r""" Notes ----- - note that .. math:: \widetilde{\pi}^P is in corrs_pressure -> from time correctors only 'u', 'dp' are needed. """ dim = problem.domain.mesh.dim vu, vp, vn, vpp1, vppp1 = var_names vdp = 'd' + vp variables = micro_problem.create_variables() to_output = variables.create_output micro_u, micro_p = variables[vu], variables[vp] micro_coor = micro_u.field.get_coor() nodes_yc = micro_problem.domain.regions['Yc'].vertices join = os.path.join format = get_print_info(problem.domain.mesh.n_el, fill='0')[1] for ii, iel in enumerate(region.cells): print('ii: %d, iel: %d' % (ii, iel)) pressure = pressures[-1][ii, 0, 0, 0] us = corrs_pressure[vu].data * pressure add_strain_rs(corrs_rs, strain, vu, dim, ii, out=us) ut = convolve_field_scalar(corrs_time_pressure[vu], pressures, ii, ts) ut += convolve_field_sym_tensor(corrs_time_rs, dstrains, vu, dim, ii, ts) u1 = us + ut u_mic = compute_u_from_macro(strain, micro_coor, ii) + u1 ps = corrs_pressure[vp].data * pressure pt = convolve_field_scalar(corrs_time_pressure[vdp], pressures, ii, ts) pt += convolve_field_sym_tensor(corrs_time_rs, dstrains, vdp, dim, ii, ts) p_hat = ps + pt # \eta_k \partial_k^x p p1 = combine_scalar_grad(corrs_permeability, p_grad, vn, ii) p_hat_e = micro_p.field.extend_dofs(p_hat[:, None], fill_value=0.0) p_mic = compute_p_from_macro(p_grad, micro_coor, ii)[:, None] \ + p_hat_e / eps0 p_mic[nodes_yc] = p1[:, None] # (y_k + \eta_k) \partial_k^x p p_aux = combine_scalar_grad(corrs_permeability, p_grad, vn, ii, shift_coors=micro_coor[nodes_yc]) meval = micro_problem.evaluate var_p = variables[vppp1] var_p.set_data(p_aux) dvel_m1 = meval('ev_diffusion_velocity.i1.Yc(m.K, %s)' % vppp1, verbose=False, mode='el_avg', **{vppp1: var_p}) var_p = variables[vpp1] var_p.set_data(p_hat) dvel_m2 = meval('ev_diffusion_velocity.i1.Ym(m.K, %s)' % vpp1, verbose=False, mode='el_avg', **{vpp1: var_p}) * eps0 out = {} out.update(to_output(u_mic, var_info={vu: (True, vu)}, extend=True)) out[vp] = Struct(name='output_data', mode='vertex', data=p_mic, var_name=vp, dofs=micro_p.dofs) aux = extend_cell_data(dvel_m1, micro_problem.domain, 'Yc') out['dvel_m1'] = Struct(name='output_data', mode='cell', data=aux, dofs=None) aux = extend_cell_data(dvel_m2, micro_problem.domain, 'Ym') out['dvel_m2'] = Struct(name='output_data', mode='cell', data=aux, dofs=None) suffix = get_output_suffix(iel, ts, naming_scheme, format, micro_problem.output_format) micro_name = micro_problem.get_output_name(extra=suffix) filename = join(problem.output_dir, 'recovered_' + os.path.basename(micro_name)) micro_problem.save_state(filename, out=out, ts=ts)
[docs] def recover_paraflow(problem, micro_problem, region, ts, strain, dstrains, pressures1, pressures2, corrs_rs, corrs_time_rs, corrs_alpha1, corrs_time_alpha1, corrs_alpha2, corrs_time_alpha2, var_names, naming_scheme='step_iel'): dim = problem.domain.mesh.dim vu, vp = var_names vdp = 'd' + vp micro_u = micro_problem.variables[vu] micro_coor = micro_u.field.get_coor() micro_p = micro_problem.variables[vp] nodes_y1 = micro_problem.domain.regions['Y1'].vertices nodes_y2 = micro_problem.domain.regions['Y2'].vertices to_output = micro_problem.variables.create_output join = os.path.join format = get_print_info(problem.domain.mesh.n_el, fill='0')[1] for ii, iel in enumerate(region.cells): print('ii: %d, iel: %d' % (ii, iel)) p1, p2 = pressures1[-1][ii, 0, 0, 0], pressures2[-1][ii, 0, 0, 0] us = corrs_alpha1[vu].data * p1 + corrs_alpha2[vu].data * p2 add_strain_rs(corrs_rs, strain, vu, dim, ii, out=us) ut = convolve_field_scalar(corrs_time_alpha1[vu], pressures1, ii, ts) ut += convolve_field_scalar(corrs_time_alpha2[vu], pressures2, ii, ts) ut += convolve_field_sym_tensor(corrs_time_rs, dstrains, vu, dim, ii, ts) u_corr = us + ut u_mic = compute_u_from_macro(strain, micro_coor, ii) + u_corr ps = corrs_alpha1[vp].data * p1 + corrs_alpha2[vp].data * p2 pt = convolve_field_scalar(corrs_time_alpha1[vdp], pressures1, ii, ts) pt += convolve_field_scalar(corrs_time_alpha2[vdp], pressures2, ii, ts) pt += convolve_field_sym_tensor(corrs_time_rs, dstrains, vdp, dim, ii, ts) p_corr = ps + pt p_mic = micro_p.field.extend_dofs(p_corr[:, nm.newaxis]) p_mic[nodes_y1] = p1 p_mic[nodes_y2] = p2 out = {} out.update(to_output(u_mic, var_info={vu: (True, vu)}, extend=True)) out[vp] = Struct(name='output_data', mode='vertex', data=p_mic, var_name=vp, dofs=micro_p.dofs) suffix = get_output_suffix(iel, ts, naming_scheme, format, micro_problem.output_format) micro_name = micro_problem.get_output_name(extra=suffix) filename = join(problem.output_dir, 'recovered_' + micro_name) micro_problem.save_state(filename, out=out, ts=ts)
[docs] def save_recovery_region(mac_pb, rname, filename=None): filename = get_default(filename, os.path.join(mac_pb.output_dir, 'recovery_region.vtk')) region = mac_pb.domain.regions[rname] # Save recovery region characteristic function. out = {} mask = region.get_charfun(by_cell=False, val_by_id=False) out['vmask'] = Struct(name='output_data', mode='vertex', data=mask[:, nm.newaxis], dofs=None) mask = region.get_charfun(by_cell=True, val_by_id=False) out['cmask'] = Struct(name='output_data', mode='cell', data=mask[:, nm.newaxis, nm.newaxis, nm.newaxis], dofs=None) mac_pb.save_state(filename, out=out)
_recovery_global_dict = {}
[docs] def destroy_pool(): if 'pool' in _recovery_global_dict: _recovery_global_dict['pool'].close()
atexit.register(destroy_pool) def _recovery_hook(args): idx, label, local_macro, verbose = args pb, corrs, recovery_hook = _recovery_global_dict['hook_args'] output.set_output(quiet=False) output.level = label[1] output(label[0]) output.set_output(quiet=not(verbose)) if isinstance(corrs, list): corrs = corrs[idx] return recovery_hook(pb, corrs, local_macro)
[docs] def get_recovery_points(region, eps0): rcoors = region.domain.mesh.coors[region.get_entities(0), :] rcmin = nm.min(rcoors, axis=0) rcmax = nm.max(rcoors, axis=0) nn = nm.round((rcmax - rcmin) / eps0) if nm.prod(nn) == 0: raise ValueError( 'incompatible recovery region and microstructure size!') cs = [] for ii, n in enumerate(nn): cs.append(nm.arange(n) * eps0 + rcmin[ii]) ccoors = nm.empty((int(nm.prod(nn)), nn.shape[0]), dtype=nm.float64) for ii, icoor in enumerate(nm.meshgrid(*cs, indexing='ij')): ccoors[:, ii] = icoor.flatten() return ccoors
[docs] def recover_micro_hook(micro_filename, region, macro, eps0, region_mode='el_centers', eval_mode='constant', eval_vars=None, corrs=None, recovery_file_tag='', define_args=None, output_dir=None, verbose=False): """ Parameters ---------- micro_filename : str The definition file of the microproblem. region : Region or array The macroscopic region to be recovered. If array, the centers of microscopic RVEs (Representative Volume Element) are assumed to be stored in it. If Region, the RVE centers are computed according to `region_mode`, see below. macro : dict of arrays or tuples Either macroscopic values (if array) or the tuple (mode, eval_var, nodal_values) is expected. The tuple is used to evaluate the macroscopic values in given points of RVEs (see 'eval_mode`). `mode` can be 'val', 'grad', 'div', or 'cauchy_strain'. eps0 : float The size of the microstructures (RVE). region_mode : {'el_centers', 'tiled'} If 'el_centers', the RVE centers are identical to the element centers of the macroscopic FE mesh. If 'tiled', the recovered region is tiled by rescaled RVEs. eval_mode : {'constant', 'continuous'} If 'constant', the macroscopic fields are evaluated only at the RVE centers. If 'continuous', the fields are evaluated at all points of the RVE mesh. eval_vars : list of variables The list of variables use to evaluate the macroscopic fields. corrs : dict of CorrSolution The correctors for recovery. recovery_file_tag : str The tag which is appended to the output file. define_args : dict The define arguments for the microscopic problem. output_dir : str The output directory. verbose : bool The verbose terminal output. """ import sfepy.base.multiproc_proc as multi if 'micro_problem' not in _recovery_global_dict: # Create a micro-problem instance. required, other = get_standard_keywords() required.remove('equations') conf = ProblemConf.from_file(micro_filename, required, other, verbose=False, define_args=define_args) if output_dir is not None: conf.options.output_dir = output_dir recovery_hook = conf.options.get('recovery_hook', None) pb = None if recovery_hook is not None: recovery_hook = conf.get_function(recovery_hook) if corrs is None: # Coefficients and correctors coefs_filename = conf.options.get('coefs_filename', 'coefs') coefs_filename = op.join(conf.options.get('output_dir', '.'), coefs_filename + '.h5') coefs = Coefficients.from_file_hdf5(coefs_filename) corrs = \ get_correctors_from_file_hdf5(dump_names=coefs.save_names) pb = Problem.from_conf(conf, init_equations=False, init_solvers=False) _recovery_global_dict['micro_problem'] = pb, corrs, recovery_hook else: pb, corrs, recovery_hook = _recovery_global_dict['micro_problem'] is_multiproc = pb.conf.options.get('multiprocessing', True)\ and multi.use_multiprocessing if recovery_hook is not None: if eval_mode not in ['constant', 'continuous']: raise ValueError(f"Evaluation mode '{eval_mode}' not implemented!") if isinstance(region, Region): if region_mode == 'el_centers': dim = region.domain.mesh.dim ccoors = region.cmesh.get_centroids(dim) ccoors = ccoors[region.get_cells(), :] recovery_ids = region.get_entities(-1) recovery_ids_s = [f'(element {k})' for k in recovery_ids] elif region_mode == 'tiled': ccoors = get_recovery_points(region, eps0) recovery_ids = nm.arange(ccoors.shape[0]) recovery_ids_s = [f'(point {k})' for k in recovery_ids] else: raise ValueError( f"Region mode '{region_mode}' not implemented!") else: ccoors = region recovery_ids = nm.arange(ccoors.shape[0]) recovery_ids_s = [f'(point {k})' for k in recovery_ids] nrec = ccoors.shape[0] output(f'recovering {nrec} microsctructures...') timer = Timer(start=True) output_level, output_fun = output.level, output.output_function _recovery_global_dict['hook_args'] = pb, corrs, recovery_hook if is_multiproc: multi_local_macro = [] num_workers = nm.min([multi.cpu_count(), nrec]) _recovery_global_dict['pool'] = multi.Pool(processes=num_workers) else: outs = [] # Recover micro. region mesh = pb.domain.mesh bbox = mesh.get_bounding_box() mic_coors = (mesh.coors - 0.5 * (bbox[1, :] + bbox[0, :])) * eps0 coors, conn, ndoffset, rec_ids = [], [], 0, [] if eval_vars is not None and eval_mode == 'constant': new_macro = {} for k, v in macro.items(): if isinstance(v, tuple): emode, evar, val = v efield = eval_vars[evar].field new_macro[k] = efield.evaluate_at(ccoors, val, mode=emode) else: new_macro[k] = v macro = new_macro for ii in range(nrec): local_coors = mic_coors.copy() local_coors[:, :ccoors.shape[1]] += ccoors[ii] coors.append(local_coors) conn.append(mesh.get_conn(mesh.descs[0]) + ndoffset) ndoffset += mesh.n_nod rec_ids.append(nm.ones((mesh.n_el,)) * recovery_ids[ii]) label = (f'micro: {ii + 1}/{nrec} {recovery_ids_s[ii]}', output_level) local_macro = {'eps0': eps0} if eval_vars is not None and eval_mode == 'continuous': for k, v in macro.items(): if isinstance(v, tuple): emode, evar, val = v efield = eval_vars[evar].field # # Inside recovery region? # from sfepy.base.base import debug; debug() # v = nm.ones((efield.region.entities[0].shape[0], 1)) # v[evfield.vertex_remap[region.entities[0]]] = 0 # no = nm.sum(v) # aux = efield.evaluate_at(local_coors, v) # if no > 0 and (nm.sum(aux) / no) > 1e-3: # print('not inside!') # continue local_macro[k] = efield.evaluate_at(local_coors, val, mode=emode) else: local_macro[k] = v else: for k, v in macro.items(): if k.startswith('_'): local_macro[k] = v else: local_macro[k] = v[ii, ...] if is_multiproc: multi_local_macro.append((ii, label, local_macro, verbose)) else: outs.append(_recovery_hook((ii, label, local_macro, verbose))) rec_ids = nm.hstack(rec_ids).reshape((-1, 1, 1, 1)) if is_multiproc: pool = _recovery_global_dict['pool'] outs = list(pool.map(_recovery_hook, multi_local_macro)) output.level, output.output_function = output_level, output_fun # Collect output data - by region outregs_data = {} outregs_info = {None: ('ALL', nm.arange(pb.domain.mesh.n_el))} vn_reg = {None: None} for k, v in outs[0].items(): if hasattr(v, 'region_name'): rn = v.region_name if rn not in outregs_info: reg = pb.domain.regions[rn] outregs_info[rn] = (rn, reg.get_entities(-1)) else: vn = getattr(v, 'var_name', None) if vn not in vn_reg: reg = pb.create_variables(vn)[vn].field.region rn = reg.name vn_reg[vn] = rn if rn not in outregs_info: outregs_info[rn] = (rn, reg.get_entities(-1)) else: rn = vn_reg[vn] if rn in outregs_data: outregs_data[rn].append(k) else: outregs_data[rn] = [k] nrve = len(coors) coors = nm.vstack(coors) ngroups = nm.tile(mesh.cmesh.vertex_groups.squeeze(), (nrve,)) conn = nm.vstack(conn) cgroups = nm.tile(mesh.cmesh.cell_groups.squeeze(), (nrve,)) # Get region mesh and data output_dir = pb.conf.options.get('output_dir', '.') for rn in outregs_data.keys(): rlabel, cidxs = outregs_info[rn] gcidxs = nm.hstack([cidxs + mesh.n_el * ii for ii in range(nrve)]) rconn = conn[gcidxs] remap = -nm.ones((coors.shape[0],), dtype=nm.int32) remap[rconn] = 1 vidxs = nm.where(remap > 0)[0] remap[vidxs] = nm.arange(len(vidxs)) rconn = remap[rconn] rcoors = coors[vidxs, :] out = {} for ivar in outregs_data[rn]: data = [outs[ii][ivar].data for ii in range(nrve)] out[ivar] = Struct(name='output_data', mode=outs[0][ivar].mode, data=nm.vstack(data)) out['recovery_id'] = Struct(name='output_data', mode='cell', data=rec_ids[gcidxs, ...]) micro_name = pb.get_output_name(extra='recovered%s_%s' % (recovery_file_tag, rlabel)) filename = op.join(output_dir, op.basename(micro_name)) mesh_out = Mesh.from_data('recovery_%s' % rlabel, rcoors, ngroups[vidxs], [rconn], [cgroups[gcidxs]], [mesh.descs[0]]) mesh_out.write(filename, io='auto', out=out) output(f'output saved to "{filename}"') output('...done in %.2f s' % timer.stop())