Merge branch 'master' of github.com:LightForm-group/damask-parse

damask_parse/readers.py

@@ -184,12 +184,11 @@ def read_geom(geom_path):
     return geometry
 
 
-def read_spectral_stdout(path):
-
+def read_spectral_stdout(path, encoding="utf8"):
     path = Path(path)
     inc_split_str = r'\s#{75}'
 
-    with path.open('r', encoding='utf8') as handle:
+    with path.open("r", encoding=encoding) as handle:
         lines = handle.read()
 
         inc_split = re.split(inc_split_str, lines)
@@ -412,7 +411,12 @@ def read_HDF5_file(
 
     """
     if not geom_path:
-        geom_path = Path(hdf5_path).parent / 'geom.vtr'
+        sim_dir = Path(hdf5_path).parent
+        geom_path = sim_dir / 'geom.vtr' # known: v3 alpha 3
+        if not geom_path.is_file():
+            geom_path = sim_dir / "geom.vti" # known: v3 alpha 7
+        if not geom_path.is_file():
+            raise ValueError(f"Cannot find geometry file in path {sim_dir!r}")
 
     # Open DAMASK output file if required
     if operations or volume_data or phase_data or field_data or grain_data:

damask_parse/utils.py

@@ -1,5 +1,7 @@
 """`damask_parse.utils.py`"""
 
+from contextlib import contextmanager
+import os
 from pathlib import Path
 from subprocess import run, PIPE
 import copy
@@ -234,6 +236,17 @@ def format_1D_masked_array(arr, fill_symbol='x'):
     return [x.item() if not isinstance(x, np.ma.core.MaskedConstant)
             else fill_symbol for x in arr]
 
+def format_2D_masked_array(arr, fill_symbol="x"):
+    out = []
+    for x in arr:
+        sub = []
+        for i in x:
+            val = (
+                i.item() if not isinstance(i, np.ma.core.MaskedConstant) else fill_symbol
+            )
+            sub.append(val)
+        out.append(sub)
+    return out
 
 def masked_array_from_list(arr, fill_value='x'):
     """Generate a (masked) array from a 1D list whose elements may contain a fill value."""
@@ -710,7 +723,12 @@ def increment_generator(increments, sim_data):
     inc_prefix = 'increment_'
 
     for inc in parse_inc_specs_using_result_obj(increments, sim_data):
-        sim_data = sim_data.view('increments', f"{inc_prefix}{inc}")
+        try:
+            # known: v3 alpha 3
+            sim_data = sim_data.view('increments', f"{inc_prefix}{inc}")
+        except TypeError:
+            # known: v3 alpha 7
+            sim_data = sim_data.view(increments=f"{inc_prefix}{inc}")
 
         yield inc, sim_data
 
@@ -791,7 +809,12 @@ def get_phase_data(sim_data, field_name, phase_name, increments,
     phase_data = []
     incs_valid = []
     for inc, sim_data in increment_generator(increments, sim_data):
-        data = sim_data.view('phases', phase_name).get(output=field_name)
+        try:
+            # known: v3 alpha 3
+            data = sim_data.view('phases', phase_name).get(output=field_name)
+        except TypeError:
+            # known: v3 alpha 7
+            data = sim_data.view(phases=phase_name).get(output=field_name)
 
         if data is None:
             print(f"Could not find field '{field_name}' for phase "
@@ -1521,24 +1544,38 @@ def get_volume_element_materials(volume_element, homog_schemes=None, phases=None
                     msg = 'Orientation `unit_cell_alignment` must be specified.'
                     raise ValueError(msg)
 
-                if volume_element['orientations']['unit_cell_alignment'].get('y') == 'b':
-                    # Convert from y//b to x//a:
+                transform_angle = None
+                x_align = volume_element['orientations']['unit_cell_alignment'].get('x')
+                y_align = volume_element['orientations']['unit_cell_alignment'].get('y')
+                if x_align == "a*" or y_align == "b":
+                    # rotate by 30 degrees about `z`
+                    transform_angle = np.pi/6
+                elif x_align == "b*" or y_align == "a":
+                    # rotate by 90 degrees about `-z`
+                    transform_angle = -np.pi/2
+                elif x_align == "a" or y_align == "b*":
+                    # already in DAMASK convention
+                    pass 
+                elif x_align == "b" or y_align == "a*":
+                    # rotate by 120 degrees about `z`
+                    transform_angle = 2*np.pi/3
+                else:
+                    msg = (f'Cannot convert from the following specified unit cell '
+                           f'alignment to DAMASK-compatible unit cell alignment (x//a): '
+                           f'{volume_element["orientations"]["unit_cell_alignment"]}')
+                    NotImplementedError(msg)
+
+                if transform_angle is not None:
                     hex_transform_quat = axang2quat(
                         volume_element['orientations']['P'] * np.array(
                             [0, 0, 1], dtype=np.longdouble),
-                        np.longdouble(np.pi/6)
+                        np.longdouble(transform_angle)
                     )
                     mat_i_const_j_ori = multiply_quaternions(
                         q1=hex_transform_quat,
                         q2=mat_i_const_j_ori,
                         P=volume_element['orientations']['P'],
-                    )
-
-                elif volume_element['orientations']['unit_cell_alignment'].get('x') != 'a':
-                    msg = (f'Cannot convert from the following specified unit cell '
-                           f'alignment to DAMASK-compatible unit cell alignment (x//a): '
-                           f'{volume_element["orientations"]["unit_cell_alignment"]}')
-                    NotImplementedError(msg)
+                    )  
 
             if volume_element['orientations']['P'] != P:
                 # Make output quaternions consistent with desired "P" convention, as
@@ -1790,3 +1827,105 @@ def spread_orientations(volume_element, phase_names, sigmas):
                 ]
 
     return volume_element
+
+@contextmanager
+def working_directory(path):
+    """Change to a working directory and return to previous working directory on exit."""
+    prev_cwd = Path.cwd()
+    os.chdir(path)
+    try:
+        yield
+    finally:
+        os.chdir(prev_cwd)
+
+def visualise_static_outpurts(outputs, result, parsed_outs):
+    """Create separate VTK file for grain and phase maps."""
+
+    static_outputs = ['grain', 'phase']
+
+    if isinstance(outputs, list):
+        static_outputs = list(set(outputs).intersection(static_outputs))
+        if len(static_outputs) > 0:
+            v = result.geometry0
+
+            for static_output in static_outputs:
+                outputs.remove(static_output)
+                dat_array = parsed_outs['field_data'][static_output]['data']
+                try:
+                    # known: v3 alpha 3
+                    v.add(dat_array.flatten(order='F'), label=static_output)
+                except AttributeError:
+                    # known: v3 alpha 7
+                    v.set(data=dat_array.flatten(order='F'), label=static_output)
+
+            v.save('static_outputs')
+
+    return outputs
+
+def generate_viz(hdf5_path, viz_spec, parsed_outs):
+    if viz_spec is not None:
+
+        from damask import Result
+
+        if viz_spec is True:
+            viz_spec = [{}]
+        elif isinstance(viz_spec, dict):
+            viz_spec = [viz_spec]
+
+        os.mkdir('viz')
+        with working_directory('viz'):
+
+            result = Result(hdf5_path)
+
+            for viz_dict_idx, viz_dict in enumerate(viz_spec, 1):
+
+                if len(viz_spec) > 1:
+                    viz_dir = str(viz_dict_idx)
+                    os.mkdir(viz_dir)
+                else:
+                    viz_dir = '.'
+                with working_directory(viz_dir):
+
+                    try:
+                        # all incs if not specified:
+                        incs_spec = viz_dict.get('increments', None)
+                        parsed_incs = parse_inc_specs_using_result_obj(incs_spec, result)
+                        try:
+                            # known: v3 alpha 3
+                            result = result.view('increments', parsed_incs)
+                        except TypeError:
+                            # known: v3 alpha 7
+                            result = result.view(increments=parsed_incs)
+
+                        # all phases if not specified:
+                        phases = viz_dict.get('phases', True)
+                        try:
+                            # known: v3 alpha 3
+                            result = result.view('phases', phases)
+                        except TypeError:
+                            # known: v3 alpha 7
+                            result = result.view(phases=phases)
+
+                        # all homogs if not specified:
+                        homogs = viz_dict.get('homogenizations', True)
+                        try:
+                            # known: v3 alpha 3
+                            result = result.view('homogenizations', homogs)
+                        except TypeError:
+                            # known: v3 alpha 7
+                            result = result.view(homogenizations=homogs)
+
+                        # all outputs if not specified:
+                        outputs = viz_dict.get('fields', '*')
+
+                        outputs = visualise_static_outpurts(outputs, result, parsed_outs)
+
+                        # result.save_VTK(output=outputs)
+                        result.export_VTK(output=outputs)
+
+                    except Exception as err:
+                        print(f'Could not save VTK files for visualise item: {viz_dict}. '
+                              f'Exception was: {err}')
+                        continue
+
+

damask_parse/writers.py

@@ -8,6 +8,7 @@
 from ruamel.yaml import YAML
 
 from damask_parse.utils import (
+    format_2D_masked_array,
     zeropad,
     format_1D_masked_array,
     align_orientations,
@@ -78,7 +79,8 @@ def write_geom(dir_path, volume_element, name='geom.vtr'):
     return geom_path
 
 
-def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None, name='load.yaml'):
+def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None,
+                    name='load.yaml', write_2D_arrs=False):
     """Write the load file for a DAMASK simulation.
 
     Parameters
@@ -90,6 +92,10 @@ def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None,
     initial_conditions : dict of (str: dict)
     name : str, optional
         Name of the load file to write. By default, set to "load.yaml".
+    write_2D_arrs : bool
+        If True, write out mechanical boundary condition arrays (stress, deformation
+        gradient, etc.) as list of lists (known to work in v3.0.0-alpha7) rather than 1D
+        lists (known to work in v3.0.0-alpha3).
 
     Returns
     -------
@@ -105,6 +111,10 @@ def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None,
         }
 
     load_steps = []
+    if write_2D_arrs:
+        fmt_arr_func = format_2D_masked_array
+    else:
+        fmt_arr_func = lambda x: format_1D_masked_array(x.flat)
 
     for load_case in load_cases:
 
@@ -179,7 +189,7 @@ def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None,
                        'passed as a masked array.')
                 raise ValueError(msg)
 
-            bc_mech[dg_arr_sym] = format_1D_masked_array(dg_arr.flat)
+            bc_mech[dg_arr_sym] = fmt_arr_func(dg_arr)
 
         else:
             if isinstance(stress_arr, np.ma.core.MaskedArray):
@@ -203,16 +213,16 @@ def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None,
                                '`vel_grad`')
                         raise ValueError(msg)
 
-                bc_mech[dg_arr_sym] = format_1D_masked_array(dg_arr.flat)
-                bc_mech[stress_arr_sym] = format_1D_masked_array(stress_arr.flat)
+                bc_mech[dg_arr_sym] = fmt_arr_func(dg_arr)
+                bc_mech[stress_arr_sym] = fmt_arr_func(stress_arr)
 
             else:
                 if dg_arr is not None:
                     msg = ('To use mixed boundary conditions, `stress` or `stress_rate`'
                            f'must be passed as a masked array.')
                     raise ValueError(msg)
 
-                bc_mech[stress_arr_sym] = format_1D_masked_array(stress_arr.flat)
+                bc_mech[stress_arr_sym] = fmt_arr_func(stress_arr)
 
         if rot is not None:
             rot = np.array(rot)
@@ -246,7 +256,8 @@ def write_load_case(dir_path, load_cases, solver=None, initial_conditions=None,
     dir_path = Path(dir_path).resolve()
     load_path = dir_path.joinpath(name)
     yaml = YAML()
-    yaml.dump(load_data, load_path)
+    with load_path.open("wt", newline="\n") as fp:
+        yaml.dump(load_data, fp)
 
     return load_path
 
@@ -376,7 +387,8 @@ def write_material(homog_schemes, phases, volume_element, dir_path, name='materi
     dir_path = Path(dir_path).resolve()
     mat_path = dir_path.joinpath(name)
     yaml = YAML()
-    yaml.dump(mat_data_fmt, mat_path)
+    with mat_path.open("wt", newline="\n") as fp:
+        yaml.dump(mat_data_fmt, fp)
 
     return mat_path