#=============================================================================
#
#  Copyright (c) Kitware, Inc.
#  All rights reserved.
#  See LICENSE.txt for details.
#
#  This software is distributed WITHOUT ANY WARRANTY; without even
#  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
#  PURPOSE.  See the above copyright notice for more information.
#
#=============================================================================
from __future__ import print_function
import datetime
import os
print('loading', os.path.basename(__file__))

import smtk
import smtk.attribute
import smtk.model
from cardformat import CardFormat
from namelist import Namelist

# ---------------------------------------------------------------------
# Strings used in "if_condition" arguments (for convenience only)
NAMELIST = 'namelist'
THERMAL_ANALYSIS = 'thermal-analysis'
ONLY_THERMAL_ANALYSIS = 'only-thermal-analysis'
FLOW_ANALYSIS = 'flow-analysis'
INVISCID_FLOW = 'inviscid_flow'
VISCOUS_FLOW = 'viscous-flow'
FLUID_PHASE = 'fluid-phase'
SOLID_PHASE = 'solid-phase'
MASS_LIMITER = 'mass-limiter'
BC_INFLOW = 'bc-inflow'
VOID_MATERIAL = 'void'  # (must match material/material-type item)
ENCLOSURE_RADIATION = 'enclosure-radiation'
MOVING_RADIATION = 'moving-radiation'
INDUCTION_HEATING = 'induction heating'

# ---------------------------------------------------------------------
class Writer:
  '''Top level writer class for Truchas input files.
  '''

# ---------------------------------------------------------------------
  def __init__(self, operator_spec, logger,
    mesh_filename='NOT-FOUND',
    altmesh_filename='NOT-FOUND'):
    '''
    '''
    self.sim_atts = smtk.attribute.Resource.CastTo(operator_spec.find('attributes').value())
    print('sim_atts', self.sim_atts)
    if self.sim_atts is None:
        msg = 'ERROR - No simulation attributes'
        print(msg)
        raise Exception(msg)

    model_entity = smtk.model.Entity.CastTo(operator_spec.find('model').objectValue(0))
    self.model_resource = smtk.model.Resource.CastTo(model_entity.resource())

    if self.model_resource is None:
        msg = 'ERROR - No model'
        print(msg)
        raise Exception(msg)

    self.analyses = list()
    self.categories = list()
    self.background_material_id = 0
    self.format_table = dict()
    self.interface_set_ids = list()
    self.enclosure_surface_set_ids = list()
    self.moving_enclosure_surface_set_ids = list()
    self.material_number_dict = {0: None}  # <phase att id, material number>
    self.mesh_file = mesh_filename
    self.altmesh_file = altmesh_filename
    self.namelist_sequence = list()
    self.skip_void_material = False
    self.out = None

# ---------------------------------------------------------------------
  def write(self, path, namelist_sequence, format_table):
    '''
    '''
    self.namelist_sequence = namelist_sequence
    self.format_table = format_table

    self._setup()

    completed = False
    with open(path, 'w') as out:
      # Write any user-provided description
      desc_att = self.sim_atts.findAttribute('description')
      if desc_att:
        desc = desc_att.findString('description').value(0)
        if desc:
          out.write('{}\n'.format(desc))
        else:
          out.write('* Truchas input file\n')

      # Add cmb output strings
      dt_string = datetime.datetime.now().strftime('%d-%b-%Y  %H:%M')
      out.write('* Generated by CMB {}\n'.format(dt_string))

      analysis_string = ', '.join(self.analyses)
      out.write('* Analysis Modules: {}\n'.format(analysis_string))


      self.out = out

      for namelist in self.namelist_sequence:
        # Check conditions
        if not CardFormat.test_conditions(namelist.if_condition):
          continue

        if namelist.separator:
          out.write('\n')
          label = '### {} '.format(namelist.separator)
          line = label.ljust(80, '#')
          out.write(line)
          out.write('\n')

        format_list = self.format_table.get(namelist.title)

        # Namelists can assign custom method
        if namelist.custom_method is not None:
          if not hasattr(self, namelist.custom_method):
            print('ERROR: For namelist', namelist.title, \
              ', custom_method', namelist.custom_method, \
              'not found')
          else:
            method = getattr(self, namelist.custom_method)
            method(namelist, format_list)
          continue

        # Else use the default namelist writer
        else:
          self._write_namelist_default(namelist, format_list)

      completed = True
      print('Wrote project file %s' % path)
    return completed

# ---------------------------------------------------------------------
  def _start_namelist(self, namelist_or_title):
    if isinstance(namelist_or_title, Namelist):
      title = namelist_or_title.title
    else:
      title = namelist_or_title
    self.out.write('\n&%s\n' % title)

# ---------------------------------------------------------------------
  def _finish_namelist(self):
    self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_namelist_default(self, namelist, format_list):
    '''
    '''
    print('Writing', namelist.title, 'namelists')

    # If namelist specifies attribute, process each one
    if namelist.att_type is not None:
      att_list = self.sim_atts.findAttributes(namelist.att_type)
      att_list.sort(key=lambda att: att.name())
      for att in att_list:
        if att.isMemberOf(self.categories):
          self._write_default_att(att, namelist, format_list)
      return

    # Otherwise process each card in format_list
    else:
      self.out.write('\n&%s\n' % namelist.title)
      for card in format_list:
        if card.att_type is None:
          raise RuntimeError('card {} missing att_type'.format(card.keyword))
        att_list = self.sim_atts.findAttributes(card.att_type)
        for att in att_list:
          if att.isMemberOf(self.categories):
            card.write(self.out, att, base_item_path=namelist.base_item_path)
      self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_default_att(self, att, namelist, format_list):
    '''Writes namelist for 1 attribute
    '''
    self.out.write('\n&%s\n' % namelist.title)

    for card in format_list:
      if card.att_type is None:
        card.write(
          self.out, att, base_item_path=namelist.base_item_path)
      else:
        card_att_list = self.sim_atts.findAttributes(card.att_type)
        for card_att in card_att_list:
          card.write(self.out, card_att)
    self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_mesh(self, namelist, format_list):
    '''
    '''
    print('Writing', namelist.title)
    self._start_namelist(namelist)

    att_list = self.sim_atts.findAttributes(namelist.att_type)
    mesh_att = att_list[0]
    mesh_comment = None
    for card in format_list:
      if not card.is_custom:
        card.write(self.out, mesh_att)
      elif 'mesh_file' == card.keyword:
        mesh_filename = os.path.basename(self.mesh_file)
        CardFormat.write_value(self.out, card.keyword, mesh_filename)
        mesh_comment = self.mesh_file
      elif 'altmesh_file' == card.keyword:
        mesh_filename = os.path.basename(self.altmesh_file)
        CardFormat.write_value(self.out, card.keyword, mesh_filename)
        mesh_comment = self.altmesh_file
      elif 'interface_side_sets' == card.keyword:
        if self.interface_set_ids:
          side_set_string = ','.join([str(x) for x in self.interface_set_ids])
          CardFormat.write_value(
            self.out, 'interface_side_sets', side_set_string, quote_string=False)

    self._finish_namelist()

    if mesh_comment:
      self.out.write('* Mesh file on local file system at:\n')
      self.out.write('*  {}\n'.format(mesh_comment))

# ---------------------------------------------------------------------
  def _write_physics(self, namelist, format_list):
    """Need custom logic because analysis attribute has no categories

    """
    print('Writing', namelist.title)
    self.out.write('\n&%s\n' % namelist.title)
    for card in format_list:
      card.write(self.out, None)
    self.out.write('/\n')


# ---------------------------------------------------------------------
  def _write_enclosure(self, namelist, format_list):
    '''
    '''
    if not 'Enclosure Radiation' in self.categories:
      return
    print('Writing', namelist.title)
    self.out.write('\n&%s\n' % namelist.title)

    att_list = self.sim_atts.findAttributes(namelist.att_type)
    enclosure_att = att_list[0]
    for card in format_list:
      if not card.is_custom:
        if card.att_type is None:
          card.write(self.out, enclosure_att, base_item_path=namelist.base_item_path)
        else:
            card_att_list = self.sim_atts.findAttributes(card.att_type)
            card_att = card_att_list[0]
            card.write(self.out, card_att)
      elif 'mesh_file' == card.keyword:
        CardFormat.write_value(self.out, card.keyword, self.mesh_file)
      elif 'side_set_ids' == card.keyword:
        side_set_string = ', '.join([str(x) for x in self.enclosure_surface_set_ids])
        CardFormat.write_value(
          self.out, card.keyword, side_set_string, quote_string=False)
      elif 'displacement_set_ids' == card.keyword:
        side_set_string = ', '.join([str(x) for x in self.moving_enclosure_surface_set_ids])
        CardFormat.write_value(
          self.out, card.keyword, side_set_string, quote_string=False)
      elif 'symmetries' == card.keyword:
        full_path = namelist.base_item_path + '/' + card.item_path
        group_item = enclosure_att.itemAtPath(full_path, '/')
        number_groups = group_item.numberOfGroups()
        if number_groups < 1:
          continue

        symmetry_list = list()
        for i in range(number_groups):
          string_item = group_item.item(i, 0)
          symmetry_value = string_item.value(0)
          if symmetry_value.startswith('Rot'):
            folds_item = string_item.activeChildItem(0)
            folds_value = folds_item.value(0)
            symmetry_value = '%s%d' % (symmetry_value, folds_value)
          symmetry_list.append(symmetry_value)
        combined_value = ", ".join(symmetry_list)
        CardFormat.write_value(
          self.out, card.keyword, combined_value, quote_string=False)

      elif 'displacement_sequence' == card.keyword:
        full_path = namelist.base_item_path + '/' + card.item_path
        group_item = enclosure_att.itemAtPath(full_path, '/')
        CardFormat.write_value(
          self.out, card.keyword, '', quote_string=False, tab=13)

        number_groups = group_item.numberOfGroups()
        last_index = number_groups - 1
        for i in range(number_groups):
          double_item = group_item.item(i, 0)
          x = double_item.value(1)
          y = double_item.value(2)
          z = double_item.value(3)
          line = '    %f, %f, %f' % (x, y, z)
          if i != last_index:
            line += ','
          line += '\n'
          self.out.write(line)
      else:
        print('WARNING skipping custom card for', card.keyword)

    self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_enclosure_radiation(self, namelist, format_list):
    '''
    '''
    if not 'Enclosure Radiation' in self.categories:
      return
    print('Writing', namelist.title)
    self.out.write('\n&%s\n' % namelist.title)

    att_list = self.sim_atts.findAttributes(namelist.att_type)
    enclosure_att = att_list[0]
    for card in format_list:
      if not card.is_custom:
        card.write(self.out, enclosure_att, base_item_path=namelist.base_item_path)
      elif 'time_sequence' == card.keyword:
        full_path = namelist.base_item_path + '/' + card.item_path
        group_item = enclosure_att.itemAtPath(full_path, '/')
        CardFormat.write_value(
          self.out, card.keyword, '', quote_string=False, tab=13)

        number_groups = group_item.numberOfGroups()
        last_index = number_groups - 1
        for i in range(number_groups):
          double_item = group_item.item(i, 0)
          time_value = double_item.value(0)
          line = '    %f' % time_value
          if i != last_index:
            line += ','
          line += '\n'
          self.out.write(line)
      else:
        print('WARNING skipping custom card for', card.keyword)

    self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_thermal_condition(self, namelist, format_list):
    '''Common writer for ds_boundary_condition, ds_interface_condition, enclosure-surface
    '''
    if not 'Heat Transfer' in self.categories:
      return
    print('Writing %s namelists' % (namelist.title))

    att_list = self.sim_atts.findAttributes(namelist.att_type)
    att_list.sort(key=lambda att: att.name())
    for att in att_list:
      #print('Writing att ', att.name())
      reference_item = att.associations()
      if reference_item is None or (0 == reference_item.numberOfValues()):
        print('Skipping attribute type \"%s\", name \"%s\" -- no associations' % \
          (att.type(), att.name()))
        continue

      self._start_namelist(namelist)
      CardFormat.write_value(self.out, 'name', att.name())
      CardFormat.write_value(self.out, 'variable', 'temperature')

      # To get condition, split the type
      split = att.type().split('.')
      CardFormat.write_value(self.out, 'condition', split[-1])

      # Write data
      data_constant_list = list()
      for i in range(att.numberOfItems()):
        item = att.item(i)
        if item.isExpression():
          expression_att = item.expression(0)
          CardFormat.write_value(self.out, 'data_function', expression_att.name())
        else:
          data_constant_list.append(item.value(0))
      if data_constant_list:
        string_list = ['{:}'.format(val) for val in data_constant_list]
        data_string = ', '.join(string_list)
        CardFormat.write_value(
          self.out, 'data_constant', data_string, quote_string=False)

      face_set_string = CardFormat.get_model_entity_ids(reference_item, as_string=True)
      CardFormat.write_value(
        self.out, 'face_set_ids', face_set_string, quote_string=False)

      self._finish_namelist()
    return

      # # Logic for enclosure_surface is straightforward
      # if 'enclosure-surface' == condition_type:
      #   print('Call write_default')
      #   self._write_default_att(att, namelist, format_list)
      #   continue

      # # Common logic for ds-boundary-condition and ds-interface-condition
      # condition_item = att.itemAtPath(namelist.base_item_path, '/')
      # #print('  condition_item', condition_item, condition)item.name())
      # condition_value = condition_item.value(0)
      # num_conditions = condition_item.numberOfActiveChildrenItems()
      # for i in range(num_conditions):
      #   active_child = condition_item.activeChildItem(i)
      #   if not active_child.isEnabled():
      #     continue

      #   active_item = active_child
      #   output_name = att.name()
      #   if num_conditions > 1:
      #     # For htc-plus-radition, override the condition value
      #     condition_value = active_item.name()
      #     # And create separate BC names
      #     output_name += '::%s' % condition_value

      #   self.out.write('\n&%s\n' % namelist.title)
      #   CardFormat.write_value(self.out, 'name', output_name)

      #   # Write the condition card
      #   CardFormat.write_value(self.out, 'condition', condition_value)

      #   # Write face set ids
      #   face_set_string = CardFormat.get_model_entity_ids(
      #     model_ent_item, as_string=True)
      #   CardFormat.write_value(
      #     self.out, 'face_set_ids', face_set_string, quote_string=False)

      #   if active_item.type() == smtk.attribute.Item.DoubleType:
      #     CardFormat.write_value(self.out, 'data_constant', active_item.value(0))
      #   if active_item.type() == smtk.attribute.Item.GroupType:
      #     value_list = list()
      #     for i in range(active_item.numberOfItemsPerGroup()):
      #       value_item = active_item.item(0, i)
      #       value_list.append(value_item.value(0))

      #     string_list = [str(x) for x in value_list]
      #     string_value = ', '.join(string_list)
      #     CardFormat.write_value(
      #       self.out, 'data_constant', string_value, quote_string=False)

      #   self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_bc(self, namelist, format_list):
    '''Common writer for fluid boundary conditions

    '''
    print('Writing', namelist.title, 'namelists')
    att_list = self.sim_atts.findAttributes(namelist.att_type)
    att_list.sort(key=lambda att: att.name())
    for att in att_list:
      if not att.isMemberOf(self.categories):
        continue

      reference_item = att.associations()
      if reference_item is None or (0 == reference_item.numberOfValues()):
        print('Skipping attribute type \"%s\", name \"%s\" -- no associations' % \
          (att.type(), att.name()))
        continue

      bc_value = None
      bc_expression = None
      # Get bc_type
      bc_type = 'dirichlet'  # default
      var_item = att.findString('variable')
      if ('velocity' == var_item.value(0)):
        #item = var_item.findChild('velocity-bc-type', smtk.attribute.SearchStyle.ACTIVE_CHILDREN)
        bc_type_item = att.itemAtPath('variable/velocity-bc-type', '/')
        bc_type = bc_type_item.value(0)

        if bc_type:
          bc_value_item = att.itemAtPath(
            'variable/velocity-bc-type/velocity-group/velocity-value', '/')
          value_list = [
            bc_value_item.value(1),
            bc_value_item.value(2),
            bc_value_item.value(3)]
          string_list = [str(x) for x in value_list]
          bc_value = ', '.join(string_list)
      else:
        bc_value_item = att.itemAtPath('variable/pressure-value', '/')
        if bc_value_item.isExpression(0):
          expression_att = bc_value_item.expression(0)
          bc_expression = expression_att.name()
        else:
          bc_value = bc_value_item.value(0)

      # Get enabled state for inflow group
      inflow_item = att.itemAtPath('variable/inflow', '/')
      if inflow_item.isEnabled():
        CardFormat.Conditions.add(BC_INFLOW)

      self.out.write('\n&%s\n' % namelist.title)
      for card in format_list:
        if 'bc_type' == card.keyword:
          CardFormat.write_value(self.out, card.keyword, bc_type)
        elif 'bc_value' == card.keyword:
          if bc_expression is not None:
            CardFormat.write_value(
              self.out, card.expression_keyword, bc_expression)
          elif bc_value is not None:
            CardFormat.write_value(
              self.out, card.keyword, bc_value, quote_string=False)
        else:
          card.write(
            self.out, att, base_item_path=namelist.base_item_path)
      self.out.write('/\n')
      CardFormat.Conditions.discard(BC_INFLOW)

# ---------------------------------------------------------------------
  def _write_body(self, namelist, format_list):
    '''Common writer for BODY namelists

    '''
    print('Writing', namelist.title, 'namelists')

    # FYI the attribute type for writing body namelists is "phase"
    att_list = self.sim_atts.findAttributes(namelist.att_type)
    if not att_list:
      return

    # Local method to find attribute of a specified type on
    # a given model entity
    def _get_attribute(att_resource, model_entity, att_type):
      # Get set of attributes associated to model_entity
      ent_set = att_resource.attributes(model_entity)
      # Get list of attributes of given type
      type_list = att_resource.findAttributes(att_type)
      # Intersection gives set that meet both cases
      common_set = ent_set.intersection(set(type_list))
      if common_set:
        return common_set.pop()  # only return 1 att
      # (else)
      return None

    # Traverse phase attributes
    att_list.sort(key=lambda att: att.name())
    for att in att_list:
      if not att.isMemberOf(self.categories):
        continue

      reference_item = att.associations()
      if reference_item is None or (0 == reference_item.numberOfValues()):
        print('Skipping attribute type \"%s\", name \"%s\" -- no associations' % \
          (att.type(), att.name()))
        continue

      # Write separate BODY namelist for each model item
      for i in range(reference_item.numberOfValues()):
        model_ent = reference_item.objectValue(i)

        self._start_namelist(namelist)
        for card in format_list:
          if 'material_number' == card.keyword:
            material_number = self.material_number_dict.get(att.id())
            CardFormat.write_value(self.out, card.keyword, material_number)
          elif card.keyword == 'mesh_material_number':
            model_ent_id = CardFormat.get_model_entity_id(model_ent)
            CardFormat.write_value(self.out, card.keyword, model_ent_id)
          elif card.keyword in ['temperature', 'velocity']:
            # Initial conditions are separate attributes on the same model entity
            ic_att = _get_attribute(self.sim_atts, model_ent, card.att_type)
            if ic_att:
              card.write(self.out, ic_att)
          else:
            card.write(self.out, att)
        self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_outputs(self, namelist, format_list):
    '''Writes OUTPUTS namelist

    Needs custom logic because of item organization
      attribute type "outputs"
        double item "start-time"
        double item "end-time"
        double item  "output-dt"
        optional extensible group item "output-times"
          double item "time"
            component 0 = output_t
            component 1 = output_dt

    '''
    print('Writing', namelist.title)
    att_list = self.sim_atts.findAttributes(namelist.att_type)
    time_list = list()
    dt_list = list()
    try:
      att = att_list[0]
      start_time = att.findDouble('start-time').value(0)
      end_time = att.findDouble('end-time').value(0)
      dt = att.findDouble('output-dt').value(0)

      time_list.append(start_time)
      dt_list.append(dt)

      # Check for optional output-times
      group_item = att.itemAtPath('output-times', '/')
      num_groups = group_item.numberOfGroups()
      for i in range(num_groups):
        double_item = group_item.item(i, 0)
        time_list.append(double_item.value(0))
        dt_list.append(double_item.value(1))
      time_list.append(end_time)
    except Exception as ex:
      print('ERROR ', ex)
      return

    # Write each component as separate array
    self.out.write('\n&%s\n' % namelist.title)

    formatted_list = ['%s' % value for value in time_list]
    formatted_string = ', '.join(formatted_list)
    CardFormat.write_value(
      self.out, 'output_t', formatted_string, quote_string=False, tab=10)

    formatted_list = ['%s' % value for value in dt_list]
    formatted_string = ', '.join(formatted_list)
    CardFormat.write_value(
      self.out, 'output_dt', formatted_string, quote_string=False, tab=10)

    self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_simcontrol(self, namelist, format_list):
    '''Common writer for SIMULATION_CONTROL namelist

    '''
    print('Writing', namelist.title, 'namelist')
    att_list = self.sim_atts.findAttributes(namelist.att_type)
    att = att_list[0]

    # Get group item, which is optional
    item = att.itemAtPath('simulation-control', '/')
    if not item.isEnabled():
      return

    # (else)
    self.out.write('\n&%s\n' % namelist.title)
    for card in format_list:
      card.write(self.out, att)
    self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_function(self, namelist, format_list):
    '''Common writer for FUNCTION namelists

    '''
    print('Writing {} namelists ({})'.format(namelist.title, namelist.att_type))
    att_list = self.sim_atts.findAttributes(namelist.att_type)
    att_list.sort(key=lambda att: att.name())
    for att in att_list:
      self.out.write('\n&%s\n' % namelist.title)
      for card in format_list:
        card.write(self.out, att)

      # Write type-specific function  data
      type_item = att.findString('type')
      function_type = type_item.value(0)

      if function_type == 'polynomial':
        group_item = att.itemAtPath('type/polynomial-terms')
        num_groups = group_item.numberOfGroups()

        if num_groups < 1:
          print('Warning: function {} contents are empty -- skipping'.format(att.name()))
          continue

        # Get number of independent vars
        exp_item0 = group_item.item(0, 1)
        num_vars = exp_item0.numberOfValues()

        coef_list = [0.0] * num_groups
        # Create 2D list for exponentials [group][variable]
        exp_lists = [ [0 for x in range(num_vars)] for y in range(num_groups) ]

        for i in range(num_groups):
          coef_item = group_item.item(i, 0)
          coef_list[i] = coef_item.value()

          exp_item = group_item.item(i, 1)
          for j in range(num_vars):
            exp_lists[i][j] = exp_item.value(j)

        CardFormat.write_value(self.out, 'poly_coefficients', coef_list)
        for k in range(num_groups):
          keyword = 'poly_exponents(:, {})'.format(k+1)
          CardFormat.write_value(self.out, keyword, exp_lists[k])

        # Finish with the center value
        CardFormat('poly_refvars', item_path='type/center').write(self.out, att)
      elif function_type == 'tabular':
        CardFormat.write_value(self.out, 'tabular_dim', 1)

        group_item = att.itemAtPath('type/tabular-data')
        num_terms = group_item.numberOfGroups()
        xvals = [0.0] * num_terms
        yvals = [0.0] * num_terms
        for i in range(num_terms):
          xvals[i] = group_item.item(i, 0).value()
          yvals[i] = group_item.item(i, 1).value()
        CardFormat.write_value(self.out, 'tabular_data(1,:)', xvals)
        CardFormat.write_value(self.out, 'tabular_data(2,:)', yvals)

        # Finish with interpolation & extrapolation
        CardFormat('tabular_interp', item_path='type/interpolation').write(self.out, att)
        CardFormat('tabular_extrap', item_path='type/extrapolation').write(self.out, att)

      self.out.write('/\n')

# ---------------------------------------------------------------------
  def _write_materials(self, namelist, format_list):
    '''Writes interleaved MATERIAL, PHASE, MATERIAL_SET namelists

    '''
    print('Writing MATERIAL/PHASE/MATERIAL_SYSTEM namelists')

    # First sort materials by single- vs multiple phase. Create 2 sets:
    #  single_phase_list: [phase attribute not appearing in any transition]
    #  first_transition_list: [lowest transition in multiphase material]

    # To sort materials/phases, first create these intermediate objects:
    #  lower_dict: <phase (id), transition> for phases in lower side of transition
    #  upper_dict: <phase (id), transition> for phases in upper side of transition
    # Note that we use ids for keys because python wrappers might create
    # multiple python objects for same c++ object

    def check_phase(phase_id, transition_dict, label):
      """Checks that phase not already used in transition dictionary"""
      if phase_id in transition_dict:
        phase_att = self.sim_atts.findAttribute(phase_id)
        tpl = 'ERROR: Phase attribute \"{}\" twice used as {} part of phase transitions'
        msg = tpl.format(phase_att.name(), label)
        print(msg)
        raise RuntimeError(msg)

    lower_dict = dict()  # <phrase attribute, transition attribute>
    upper_dict = dict()  # <phrase attribute, transition attribute>
    trans_atts = self.sim_atts.findAttributes('phase-transition')
    for trans_att in trans_atts:
      if not trans_att.isValid():
        tpl = 'ERROR: invalid Phase Transition attribute \"{}\"'
        msg = tpl.format(trans_att.name())
        print(msg)
        raise RuntimeError(msg)

      # Get lower & upper components and their ids
      lower_att = trans_att.findComponent('lower').value()
      upper_att = trans_att.findComponent('upper').value()

      lower_id = lower_att.id()
      upper_id = upper_att.id()

      # Check for consistency
      if lower_id == upper_id:
        tpl = 'ERROR: same phase -- \"{}\" -- used for both sides of transition \"{}\"'
        msg = tpl.format(lower.name(), att.name())
        print(msg)
        raise RuntimeError(msg)

      lower_temp = trans_att.findDouble('lower-transition-temperature').value()
      upper_temp = trans_att.findDouble('upper-transition-temperature').value()
      if upper_temp <= lower_temp:
        tpl = 'ERROR: phase transition \"{}\" upper temp <= lower temp ({} <= {})'
        msg = tpl.format(trans_att.name(), upper_temp, lower_temp)
        print(msg)
        raise RuntimeError(msg)

      check_phase(lower_id, lower_dict, 'lower')
      check_phase(upper_id, upper_dict, 'upper')

      lower_dict[lower_id] = trans_att
      upper_dict[upper_id] = trans_att

    # Create list of "materials", each item either:
    #  a phase attribute representing a single-phase material
    #  a phase-transition attribute representing the "first" transition of multi-phase material
    material_list = list()
    phase_atts = self.sim_atts.findAttributes('phase')
    for phase_att in phase_atts:
      phase_id = phase_att.id()
      in_lower = phase_id in lower_dict
      in_upper = phase_id in upper_dict
      if not in_lower and not in_upper:
        material_list.append(phase_att)  # single-phase material
      elif in_lower and not in_upper:
        material_list.append(lower_dict[phase_id])  # multi-phase material (transition)
    # print('material_list: {}'.format(material_list))

    def set_conditions(att=None):
      """Sets CardFormat conditions for given phase attribute"""
      CardFormat.Conditions.discard(VOID_MATERIAL)
      CardFormat.Conditions.discard(FLUID_PHASE)
      if att is None:
        return
      if att.type() == 'phase.void':
        CardFormat.Conditions.add(VOID_MATERIAL)
      elif att.type() == 'phase.material':
        fluid_item = att.findGroup('fluid')
        if fluid_item and fluid_item.isEnabled():
          CardFormat.Conditions.add(FLUID_PHASE)

    # Traverse material list
    # Also keep track of transitions, to make sure *all* get used
    trans_id_list = [t.id() for t in trans_atts]
    trans_id_set = set(trans_id_list)
    for att in material_list:
      if att.type() == 'phase.void' and self.skip_void_material:
          continue

      if att.type() in ['phase.material', 'phase.void']:
        # Write single phase material
        set_conditions(att)
        self._write_phase_namelist(att)
        self._write_material_namelist(att)
        self._write_material_system_namelist(att)
      elif att.type() == 'phase-transition':
        trans_att = att
        # Write multiple phase material, starting with lowest-temp phase
        lower_att = trans_att.findComponent('lower').value()
        set_conditions(lower_att)
        self._write_phase_namelist(lower_att)
        self._write_material_namelist(lower_att)
        while trans_att is not None:
          trans_id_set.discard(trans_att.id())
          upper_att = trans_att.findComponent('upper').value()
          set_conditions(upper_att)
          self._write_phase_namelist(upper_att)
          self._write_material_namelist(upper_att)

          # Get next transition
          trans_att = lower_dict.get(upper_att.id())
        self._write_material_system_namelist(att, lower_dict)

        set_conditions(None)  # resets material conditions
      else:
        raise RuntimeError('ERROR - unexpected attribute type ', att.type())

    # Any transaction atts not removed above indicate error(s)
    if trans_id_set:
      unused_atts = [self.sim_atts.findAttribute(id) for id in trans_id_set]
      unused_att_names = [att.name() for att in unused_atts]
      tpl = 'ERROR - inconsistent transitions, probably circular: {}'
      msg = tpl.format(unused_att_names)
      raise RuntimeError(msg)


# --------------------------------------------------------------------
  def _write_material_namelist(self, phase_att):
    '''Write MATERIAL namelist for phase attribute

    This method also assigns material number
    '''
    title = 'MATERIAL'
    print('Writing namelist {} for attribute {}'.format(title, phase_att.name()))
    self._start_namelist(title)
    format_list = self.format_table.get(title)

    # Assign next material number
    number = len(self.material_number_dict)
    self.material_number_dict[phase_att.id()] = number

    CardFormat.write_value(self.out, 'material_name', phase_att.name())
    CardFormat.write_value(self.out, 'material_number', number)

    if phase_att.type() == 'phase.material':
      immobile = FLUID_PHASE not in CardFormat.Conditions
      CardFormat.write_value(self.out, 'immobile', immobile, as_boolean=True)

    for card in format_list:
      card.write(self.out, phase_att)

    # Write material density (1 for material, 0 for void)
    density = 0.0 if VOID_MATERIAL in CardFormat.Conditions else 1.0
    CardFormat.write_value(self.out, 'density', density)

    # Check if this material is the background material
    if phase_att.id() == self.background_material_id:
      CardFormat.write_value(self.out, 'material_feature','background')

    self._finish_namelist()

# ---------------------------------------------------------------------
  def _write_phase_namelist(self, phase_att):
    '''Write PHASE namelist for phase attribute
    '''
    title = 'PHASE'
    print('Writing namelist {} for attribute {}'.format(title, phase_att.name()))
    self._start_namelist(title)

    name = phase_att.name()
    CardFormat.write_value(self.out, 'name', name[:31], tab=4)

    if not VOID_MATERIAL in CardFormat.Conditions:
      format_list = self.format_table.get(title)

      CardFormat.PropertyIndex = 0
      for card in format_list:
        card.write(self.out, phase_att)
    self._finish_namelist()

# ---------------------------------------------------------------------
  def _write_material_system_namelist(self, att, lower_transition_dict={}):
    '''

    att: either phase att for single-phase, or first transition att for multi-phase
    lower_transition_dict: <phase att id, transition att>
    '''
    title = 'MATERIAL_SYSTEM'
    print('Writing', title)
    self._start_namelist(title)

    # Simple case - att is single "phase" attribute
    if att.type() in ['phase.material', 'phase.void']:
      number = self.material_number_dict.get(att.id())
      name = 'material {}'.format(number)
      CardFormat.write_value(self.out, 'name', name, tab=6)
      CardFormat.write_value(self.out, 'phases', att.name(), tab=6)
      self._finish_namelist()
      return

    elif att.type() != 'phase-transition':
      raise RuntimeError('Unexpected attribute type \"{}\"'.format(att.type()))

    # Multiple phase case
    # Traverse transitions to build data lists
    phase_list = list()
    lower_temps = list()
    upper_temps = list()
    latent_heats = list()
    smoothing_radius = None
    numbers = list()  # material numbers

    trans_att = att
    lower_att = trans_att.findComponent('lower').value()
    lower_name = lower_att.name()
    number = self.material_number_dict.get(lower_att.id())
    numbers.append(number)
    phase_list.append('\"{}\"'.format(lower_name[:31]))
    while trans_att is not None:
      upper_att = trans_att.findComponent('upper').value()
      upper_name = upper_att.name()
      phase_list.append('\"{}\"'.format(upper_name[:31]))
      number = self.material_number_dict.get(upper_att.id())
      numbers.append(number)

      lower_temp_item = trans_att.findDouble('lower-transition-temperature')
      lower_temps.append(lower_temp_item.value())
      upper_temp_item = trans_att.findDouble('upper-transition-temperature')
      upper_temps.append(upper_temp_item.value())
      latent_heat = trans_att.findDouble('latent-heat')
      latent_heats.append(latent_heat.value())

      # Smoothing radius currently part of transition, but only a single value
      # applies to the material system. This code takes the smoothing radius
      # from the *first* transition that has its item enabled.
      if smoothing_radius is None:
        smoothing_radius_item = trans_att.findDouble('smoothing-radius')
        if smoothing_radius_item and smoothing_radius_item.isEnabled():
          smoothing_radius = smoothing_radius_item.value()

      # Get next transition
      trans_att = lower_transition_dict.get(upper_att.id())

    # Check that transition temps are monotonically increasing
    temp_lists = (lower_temps, upper_temps)
    temps = [val for pair in zip(*temp_lists) for val in pair]
    #print('temps', temps)
    prev_temp = temps[0]
    for i in range(1, len(temps)):
      temp = temps[i]
      if temp < prev_temp:
        tpl = 'ERROR: temperatures for material {} not monotonically increasing: {}'
        msg = tpl.format(lower_name, temps)
        print(msg)
        raise RuntimeError(msg)
      prev_temp = temp

    numbers.sort()
    number_string = '+'.join([str(n) for n in numbers])
    name = 'material {}'.format(number_string)
    CardFormat.write_value(self.out, 'name', name)
    CardFormat.write_value(self.out, 'phases', phase_list)
    CardFormat.write_value(self.out, 'transition_temps_low', lower_temps)
    CardFormat.write_value(self.out, 'transition_temps_high', upper_temps)
    CardFormat.write_value(self.out, 'latent_heat', latent_heats)
    if smoothing_radius is not None:
      CardFormat.write_value(self.out, 'smoothing_radius', smoothing_radius)

    self._finish_namelist()

# ---------------------------------------------------------------------
  def _setup(self):
    '''
    '''
    CardFormat.ModelManager = self.model_resource

    # Get analysis type(s) and set CardFormat conditions
    analysis_att = self.sim_atts.findAttribute('analysis')
    if not analysis_att:
      raise RuntimeError('Internal Error -- missing analysis att')

    # Get categories based on enabled analyses
    ER = 'Enclosure Radiation'
    FF = 'Fluid Flow'
    HT = 'Heat Transfer'
    IH = 'Induction Heating'
    SM = 'Solid Mechanics'

    self.analyses = list()
    self.categories = list()

    # Only looking at Heat Transfer and Fluid Flow initially
    er_item = None  # enclosure radition
    ff_item = analysis_att.findVoid(FF)
    ht_item = analysis_att.findGroup(HT)
    ih_item = None  # induction heating
    sm_item = analysis_att.findVoid(SM)
    analyses = self.sim_atts.analyses()
    category_set = set()
    if ht_item.isEnabled():
      ht_name = analyses.find(HT).name()
      print("*** ht_name", ht_name, "***")


      CardFormat.Conditions.add(THERMAL_ANALYSIS)
      if not ff_item.isEnabled():
        CardFormat.Conditions.add(ONLY_THERMAL_ANALYSIS)
      self.analyses.append(HT)
      category_set |= analyses.find(HT).categories()

      er_item = ht_item.find(ER)
      print('er_item:', er_item)
      if er_item.isEnabled():
        CardFormat.Conditions.add(ENCLOSURE_RADIATION)
        self.analyses.append(ER)
        category_set |= analyses.find(ER).categories()
        # Todo MOVING_RADIAITON
        raise RuntimeError('Enclosure Radiation selected but NOT yet supported')

      ih_item = ht_item.find(IH)
      if ih_item.isEnabled():
        CardFormat.Conditions.add(INDUCTION_HEATING)
        self.analyses.append(IH)
        category_set |= analyses.find(IH).categories()

    if ff_item.isEnabled():
      CardFormat.Conditions.add(FLOW_ANALYSIS)
      self.analyses.append(FF)
      category_set |= analyses.find(FF).categories(FF)
      # Todo VISCOUS_FLOW vs INVISCID_FLOW

    if sm_item.isEnabled():
      self.analyses.append(SM)
      category_set |= analyses.find(SM).categories()
      raise RuntimeError('Solid Mechanics selected but NOT yet supported')

    self.categories = list(category_set)
    print('Enabled analyses: {}'.format(self.analyses))
    print('Enabled categories: {}'.format(self.categories))
    # Iterim logic during development
    # att_list = self.sim_atts.findAttributes('ht.solver')
    # for att in att_list:
    #   analysis_item = att.findString('analysis')
    #   if analysis_item is not None:
    #     analysis = analysis_item.value(0)
    #     if analysis == 'thermal':
    #       print('Adding analysis conditions "{}" "{}"'.format(
    #         THERMAL_ANALYSIS, ONLY_THERMAL_ANALYSIS))
    #       CardFormat.Conditions.add(THERMAL_ANALYSIS)
    #       CardFormat.Conditions.add(ONLY_THERMAL_ANALYSIS)
    #     elif analysis == 'thermal-plus-fluid':
    #       print('Adding analysis conditions "{}", "{}"'.format(
    #         THERMAL_ANALYSIS, FLOW_ANALYSIS))
    #       CardFormat.Conditions.add(THERMAL_ANALYSIS)
    #       CardFormat.Conditions.add(FLOW_ANALYSIS)

    # Get the background material attribute's id
    att_list = self.sim_atts.findAttributes('background-material')
    if att_list:
      att = att_list[0]
      item = att.findComponent('background-material')
      if item.isSet(0):
        self.background_material_id = item.value(0).id()

    # Determine if void material should be written
    void_list = self.sim_atts.findAttributes('phase.void')
    if len(void_list) == 1:
      void_att = void_list[0]
      self.skip_void_material = void_att.associations().numberOfValues() == 0
    else:
      print('Warning - missing phase.void attribute')

    # Find all interface_set_ids and enclosure_surface_set_ids
    self.interface_set_ids = list()
    self.enclosure_surface_set_ids = list()
    self.moving_enclosure_surface_set_ids = list()

    att_list = self.sim_atts.findAttributes('ht.interface')
    for att in att_list:
      reference_item = att.associations()
      if reference_item:
        self.interface_set_ids += CardFormat.get_model_entity_ids(reference_item)

      # type_item = att.findString('type')
      # if type_item.value(0) == 'ds-interface-condition':
      #   self.interface_set_ids += CardFormat.get_model_entity_ids(reference_item)
      # elif type_item.value(0) == 'enclosure-surface':
      #   surface_ids = CardFormat.get_model_entity_ids(reference_item)
      #   self.enclosure_surface_set_ids += surface_ids
      #   moving_item = att.itemAtPath('type/moving', '/')
      #   if moving_item and moving_item.isEnabled():
      #     self.moving_enclosure_surface_set_ids += surface_ids

    self.interface_set_ids.sort()
    self.enclosure_surface_set_ids.sort()
    self.moving_enclosure_surface_set_ids.sort()
    print('interface_set_ids', self.interface_set_ids)
    print('enclosure_surface_set_ids', self.enclosure_surface_set_ids)
    print('moving_enclosure_surface_set_ids', self.moving_enclosure_surface_set_ids)

    # Check for moving radiation
    # att_list = self.sim_atts.findAttributes('enclosure-radiation')
    # if att_list:
    #   enclosure_att = att_list[0]
    #   CardFormat.Conditions.add(ENCLOSURE_RADIATION)

    #   moving_item = enclosure_att.itemAtPath('enclosure/moving-radiation', '/')
    #   #print('moving_item value', moving_item.value(0))
    #   if moving_item.value(0) == 1:
    #     print('Set condition', MOVING_RADIATION)
    #     CardFormat.Conditions.add(MOVING_RADIATION)
