content: #
  #           MuPIF: Multi-Physics Integration Framework
  #               Copyright (C) 2010-2015 Borek Patzak
  #
  #    Czech Technical University, Faculty of Civil Engineering,
  #  Department of Structural Mechanics, 166 29 Prague, Czech Republic
  #
  # This library is free software; you can redistribute it and/or
  # modify it under the terms of the GNU Lesser General Public
  # License as published by the Free Software Foundation; either
  # version 2.1 of the License, or (at your option) any later version.
  #
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  # Lesser General Public License for more details.
  #
  # You should have received a copy of the GNU Lesser General Public
  # License along with this library; if not, write to the Free Software
  # Foundation, Inc., 51 Franklin Street, Fifth Floor,
  # Boston, MA  02110-1301  USA
  #
  from builtins import range
  from builtins import object
  from . import Cell
  from . import FieldID
  from . import ValueType
  from . import BBox
  from . import APIError
  from . import MupifObject
  from . import Mesh
  from .Physics import PhysicalQuantities 
  from .Physics.PhysicalQuantities import PhysicalQuantity
  
  from numpy import array, arange, random, zeros
  import numpy
  import copy
  import Pyro4
  from enum import IntEnum
  import logging
  log = logging.getLogger()
  
  try:
     import cPickle as pickle  # faster serialization if available
  except:
     import pickle
  # import logging - never use it here, it causes cPickle.PicklingError: Can't pickle <type 'thread.lock'>: attribute
  # lookup thread.lock failed
  
  # debug flag
  debug = 0
  
  
  class FieldType(IntEnum):
      """
      Represent the supported values of FieldType, i.e. FT_vertexBased or FT_cellBased.
      """
      FT_vertexBased = 1
      FT_cellBased = 2
  
  
  @Pyro4.expose
  class Field(MupifObject.MupifObject, PhysicalQuantity):
      """
      Representation of field. Field is a scalar, vector, or tensorial
      quantity defined on a spatial domain. The field, however is assumed
      to be fixed at certain time. The field can be evaluated in any spatial point
      belonging to underlying domain.
  
      Derived classes will implement fields defined on common discretizations,
      like fields defined on structured/unstructured FE meshes, FD grids, etc.
  
      .. automethod:: __init__
      .. automethod:: _evaluate
      """
      def __init__(self, mesh, fieldID, valueType, units, time, values=None, fieldType=FieldType.FT_vertexBased, objectID=0, metaData={}):
          """
          Initializes the field instance.
  
          :param Mesh.Mesh mesh: Instance of a Mesh class representing the underlying discretization
          :param FieldID fieldID: Field type (displacement, strain, temperature ...)
          :param ValueType valueType: Type of field values (scalar, vector, tensor). Tensor is a tuple of 9 values. It is changed to 3x3 for VTK output automatically.
          :param Physics.PhysicalUnits units: Field value units
          :param Physics.PhysicalQuantity time: Time associated with field values
          :param values: Field values (format dependent on a particular field type, however each individual value should be stored as tuple, even scalar value)
          :type values: list of tuples representing individual values
          :param FieldType fieldType: Optional, determines field type (values specified as vertex or cell values), default is FT_vertexBased
          :param int objectID: Optional ID of problem object/subdomain to which field is related, default = 0
          :param dict metaData: Optionally pass metadata for merging
          """
          
          super(Field, self).__init__()
          self.mesh = mesh
          self.fieldID = fieldID
          self.valueType = valueType
          self.time = time
          self.uri = None  # pyro uri; used in distributed setting
          # self.log = logging.getLogger()
          self.fieldType = fieldType
          self.objectID = objectID        
          if values is None:
              if self.fieldType == FieldType.FT_vertexBased:
                  ncomponents = mesh.getNumberOfVertices()
              else:
                  ncomponents = mesh.getNumberOfCells()
              self.value = zeros((ncomponents, self.getRecordSize()))
          else:
              self.value = values
  
          if PhysicalQuantities.isPhysicalUnit(units):
              self.unit = units
          else:
              self.unit = PhysicalQuantities.findUnit(units)
              
          self.setMetadata('Units', self.unit.name())
          self.setMetadata('Type', 'mupif.Field.Field')
          self.setMetadata('Type_ID', str(self.fieldID))
          self.setMetadata('FieldType', str(fieldType))
          self.setMetadata('ValueType', str(self.valueType))
          
          self.updateMetadata(metaData)
  
      @classmethod
      def loadFromLocalFile(cls, fileName):
          """
          Alternative constructor which loads instance directly from a Pickle module.
  
          :param str fileName: File name
  
          :return: Returns Field instance
          :rtype: Field
          """
          return pickle.load(open(fileName, 'rb'))
  
      def getRecordSize(self):
          """
          Return the number of scalars per value, depending on :obj:`valueType` passed when constructing the instance.
  
          :return: number of scalars (1,3,9 respectively for scalar, vector, tensor)
          :rtype: int
          """
          if self.valueType == ValueType.Scalar:
              return 1
          elif self.valueType == ValueType.Vector:
              return 3
          elif self.valueType == ValueType.Tensor:
              return 9
          else:
              raise ValueError("Invalid value of Field.valueType (%d)." % self.valueType)
  
      def getMesh(self):
          """
          Obtain mesh.
  
          :return: Returns a mesh of underlying discretization
          :rtype: Mesh.Mesh
          """
          return self.mesh
  
      def getValueType(self):
          """
          Returns ValueType of the field, e.g. scalar, vector, tensor.
  
          :return: Returns value type of the receiver
          :rtype: ValueType
          """
          return self.valueType
  
      def getFieldID(self):
          """
          Returns FieldID, e.g. FID_Displacement, FID_Temperature.
  
          :return: Returns field ID
          :rtype: FieldID
          """
          return self.fieldID
  
      def getFieldIDName(self):
          """
          Returns name of the field.
  
          :return: Returns fieldID name
          :rtype: string
          """
          return self.fieldID.name
  
      def getFieldType(self):
          """
          Returns receiver field type (values specified as vertex or cell values)
  
          :return: Returns fieldType id
          :rtype: FieldType
          """
          return self.fieldType
  
      def getTime(self):
          """
          Get time of the field.
  
          :return: Time of field data
          :rtype: Physics.PhysicalQuantity
          """
          return self.time
  
      def evaluate(self, positions, eps=0.0):
          """
          Evaluates the receiver at given spatial position(s).
  
          :param positions: 1D/2D/3D position vectors
          :type positions: tuple, a list of tuples
          :param float eps: Optional tolerance for probing whether the point belongs to a cell (should really not be used)
          :return: field value(s)
          :rtype: Physics.PhysicalQuantity with given value or tuple of values
          """
          # test if positions is a list of positions
          if isinstance(positions, list):
              ans = []
              for pos in positions:
                  ans.append(self._evaluate(pos, eps))
              return PhysicalQuantity(ans, self.unit)
          else:
              # single position passed
              return PhysicalQuantity(self._evaluate(positions, eps), self.unit)
  
      def _evaluate(self, position, eps):
          """
          Evaluates the receiver at a single spatial position.
  
          :param tuple position: 1D/2D/3D position vector
          :param float eps: Optional tolerance
          :return: field value
          :rtype: tuple  of doubles 
  
          .. note:: This method has some issues related to https://sourceforge.net/p/mupif/tickets/22/ .
          """
          cells = self.mesh.giveCellLocalizer().giveItemsInBBox(BBox.BBox([c-eps for c in position], [c+eps for c in position]))
          # answer=None
          if len(cells):
              if self.fieldType == FieldType.FT_vertexBased:
                  for icell in cells:
                      try:
                          if icell.containsPoint(position):
                              if debug:
                                  log.debug(icell.getVertices())
                              try:
                                  answer = icell.interpolate(position, [self.value[i.number] for i in icell.getVertices()])
                              except IndexError:
                                  log.error('Field::evaluate failed, inconsistent data at cell %d' % icell.label)
                                  raise
                              return answer
  
                      except ZeroDivisionError:
                          print('ZeroDivisionError?')
                          log.debug(icell.number)
                          log.debug(position)
                          icell.debug = 1
                          log.debug(icell.containsPoint(position), icell.glob2loc(position))
  
                  log.error('Field::evaluate - no source cell found for position %s' % str(position))
                  for icell in cells:
                      log.debug(icell.number)
                      log.debug(icell.containsPoint(position))
                      log.debug(icell.glob2loc(position))
  
              else:  # if (self.fieldType == FieldType.FT_vertexBased):
                  # in case of cell based fields do compute average of cell values containing point
                  # this typically happens when point is on the shared edge or vertex
                  count = 0
                  for icell in cells:
                      if icell.containsPoint(position):
                          if debug:
                              log.debug(icell.getVertices())
  
                          try:
                              tmp = self.value[icell.number]
                              if count == 0:
                                  answer = list(tmp)
                              else:
                                  for i in answer:
                                      answer = [x+y for x in answer for y in tmp]
                              count += 1
  
                          except IndexError:
                              log.error('Field::evaluate failed, inconsistent data at cell %d' % icell.label)
                              log.error(icell.getVertices())
                              raise
                  # end loop over icells
                  if count == 0:
                      log.error('Field::evaluate - no source cell found for position %s', str(position))
                      # for icell in cells:
                      #    log.debug(icell.number, icell.containsPoint(position), icell.glob2loc(position))
                  else:
                      answer = [x/count for x in answer]
                      return answer
  
          else:
              # no source cell found
              log.error('Field::evaluate - no source cell found for position ' + str(position))
              raise ValueError('Field::evaluate - no source cell found for position ' + str(position))
  
      def getVertexValue(self, vertexID):
          """
          Returns the value associated with a given vertex.
  
          :param int vertexID: Vertex identifier
          :return: The value
          :rtype: Physics.PhysicalQuantity
          """
          if self.fieldType == FieldType.FT_vertexBased:
              return PhysicalQuantity(self.value[vertexID], self.unit)
          else:
              raise TypeError('Attempt to acces vertex value of cell based field, use evaluate instead')
          
      def getCellValue(self, cellID):
          """
          Returns the value associated with a given cell.
  
          :param int cellID: Cell identifier
          :return: The value
          :rtype: Physics.PhysicalQuantity
          """
          if self.fieldType == FieldType.FT_cellBased:
              return PhysicalQuantity(self.value[cellID], self.unit)
          else:
              raise TypeError('Attempt to acces cell value of vertex based field, use evaluate instead')
  
      def _giveValue(self, componentID):
          """
          Returns the value associated with a given component (vertex or cell).
          Depreceated, use getVertexValue() or getCellValue()
  
          :param int componentID: An identifier of a component: vertexID or cellID
          :return: The value
          :rtype: Physics.PhysicalQuantity
          """
          return PhysicalQuantity(self.value[componentID], self.unit)
      
      def giveValue(self, componentID):
          """
          Returns the value associated with a given component (vertex or cell).
  
          :param int componentID: An identifier of a component: vertexID or cellID
          :return: The value
          :rtype: tuple
          """
          return self.value[componentID]    
  
      def setValue(self, componentID, value):
          """
          Sets the value associated with a given component (vertex or cell).
  
          :param int componentID: An identifier of a component: vertexID or cellID
          :param tuple value: Value to be set for a given component, should have the same units as receiver
  
          .. Note:: If a mesh has mapping attached (a mesh view) then we have to remember value locally and record change. The source field values are updated after commit() method is invoked.
          """
          self.value[componentID] = value
  
      def commit(self):
          """
          Commits the recorded changes (via setValue method) to a primary field.
          """
      
      def getObjectID(self):
          """
          Returns field objectID.
  
          :return: Object's ID 
          :rtype: int
          """
          return self.objectID
      
      def getUnits(self):
          """
          :return: Returns units of the receiver
          :rtype: Physics.PhysicalUnits
          """
          return self.unit
  
      def merge(self, field):
          """
          Merges the receiver with given field together. Both fields should be on different parts of the domain (can also overlap), but should refer to same underlying discretization, otherwise unpredictable results can occur.
  
          :param Field field: given field to merge with.
          """
          # first merge meshes
          mesh = copy.deepcopy(self.mesh)
          mesh.merge(field.mesh)
          log.debug(mesh)
          # merge the field values
          # some type checking first
          if self.fieldType != field.fieldType:
              raise TypeError("Field::merge: fieldType of receiver and parameter is different")
          if self.fieldType == FieldType.FT_vertexBased:
              values = [0]*mesh.getNumberOfVertices()
              for v in range(self.mesh.getNumberOfVertices()):
                  values[mesh.vertexLabel2Number(self.mesh.getVertex(v).label)] = self.value[v]
              for v in range(field.mesh.getNumberOfVertices()):
                  values[mesh.vertexLabel2Number(field.mesh.getVertex(v).label)] = field.value[v]
          else:
              values = [0]*mesh.getNumberOfCells()
              for v in range(self.mesh.getNumberOfCells()):
                  values[mesh.cellLabel2Number(self.mesh.giveCell(v).label)] = self.value[v]
              for v in range(field.mesh.getNumberOfCells()):
                  values[mesh.cellLabel2Number(field.mesh.giveCell(v).label)] = field.value[v]
  
          self.mesh = mesh
          self.value = values
  
      def field2VTKData (self, name=None, lookupTable=None):
          """
          Creates VTK representation of the receiver. Useful for visualization. Requires pyvtk module.
  
          :param str name: human-readable name of the field
          :param pyvtk.LookupTable lookupTable: color lookup table
          :return: Instance of pyvtk
          :rtype: pyvtk.VtkData
          """
          import pyvtk
  
          if name is None:
              name = self.getFieldIDName()
          if lookupTable and not isinstance(lookupTable, pyvtk.LookupTable):
              log.info('ignoring lookupTable which is not a pyvtk.LookupTable instance.')
              lookupTable = None
          if lookupTable is None:
              lookupTable=pyvtk.LookupTable([(0, .231, .298, 1.0), (.4, .865, .865, 1.0), (.8, .706, .016, 1.0)], name='coolwarm')
              # Scalars use different name than 'coolwarm'. Then Paraview uses its own color mapping instead of taking
              # 'coolwarm' from *.vtk file. This prevents setting Paraview's color mapping.
              scalarsKw = dict(name=name, lookup_table='default')
          else:
              scalarsKw = dict(name=name, lookup_table=lookupTable.name)
          # see http://cens.ioc.ee/cgi-bin/cvsweb/python/pyvtk/examples/example1.py?rev=1.3 for an example
          vectorsKw = dict(name=name)  # vectors don't have a lookup_table
  
          if self.fieldType == FieldType.FT_vertexBased:
              if self.getValueType() == ValueType.Scalar:
                  return pyvtk.VtkData(self.mesh.getVTKRepresentation(), pyvtk.PointData(pyvtk.Scalars([val[0] for val in self.value], **scalarsKw), lookupTable), 'Unstructured Grid Example')
              elif self.getValueType() == ValueType.Vector:
                  return pyvtk.VtkData(self.mesh.getVTKRepresentation(), pyvtk.PointData(pyvtk.Vectors(self.value, **vectorsKw), lookupTable), 'Unstructured Grid Example')
              elif self.getValueType() == ValueType.Tensor:
                  return pyvtk.VtkData(self.mesh.getVTKRepresentation(), pyvtk.PointData(pyvtk.Tensors(self.getMartixForTensor(self.value), **vectorsKw), lookupTable), 'Unstructured Grid Example')
              
          else:
              if self.getValueType() == ValueType.Scalar:
                  return pyvtk.VtkData(self.mesh.getVTKRepresentation(), pyvtk.CellData(pyvtk.Scalars([val[0] for val in self.value], **scalarsKw), lookupTable), 'Unstructured Grid Example')
              elif self.getValueType() == ValueType.Vector:
                  return pyvtk.VtkData(self.mesh.getVTKRepresentation(), pyvtk.CellData(pyvtk.Vectors(self.value, **vectorsKw),lookupTable), 'Unstructured Grid Example')
              elif self.getValueType() == ValueType.Tensor:
                  return pyvtk.VtkData(self.mesh.getVTKRepresentation(), pyvtk.CellData(pyvtk.Tensors(self.getMartixForTensor(self.value), **vectorsKw), lookupTable), 'Unstructured Grid Example')
              
      def getMartixForTensor(self, values):
          """
          Reshape values to a list with 3x3 arrays. Usable for VTK export.
  
          :param list values: List containing tuples of 9 values, e.g. [(1,2,3,4,5,6,7,8,9), (1,2,3,4,5,6,7,8,9), ...]
          
          :return: List containing 3x3 matrices for each tensor
          :rtype: list
          """ 
          tensor = []
          for i in values:
              tensor.append(numpy.reshape(i, (3, 3)))
          return tensor
  
      def dumpToLocalFile(self, fileName, protocol=pickle.HIGHEST_PROTOCOL):
          """
          Dump Field to a file using a Pickle serialization module.
  
          :param str fileName: File name
          :param int protocol: Used protocol - 0=ASCII, 1=old binary, 2=new binary
          """
          pickle.dump(self, open(fileName, 'wb'), protocol)
  
      def field2Image2D(self, plane='xy', elevation=(-1.e-6, 1.e-6), numX=10, numY=20, interp='linear', fieldComponent=0, vertex=True, colorBar='horizontal', colorBarLegend='', barRange=(None, None), barFormatNum='%.3g', title='', xlabel='', ylabel='', fileName='', show=True, figsize=(8, 4), matPlotFig=None):
          """ 
          Plots and/or saves 2D image using a matplotlib library. Works for structured and unstructured 2D/3D fields. 2D/3D fields need to define plane. This method gives only basic viewing options, for aesthetic and more elaborated output use e.g. VTK field export with 
          postprocessors such as ParaView or Mayavi. Idea from https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html#id1
  
          :param str plane: what plane to extract from field, valid values are 'xy', 'xz', 'yz'
          :param tuple elevation: range of third coordinate. For example, in plane='xy' is grabs z coordinates in the range
          :param int numX: number of divisions on x graph axis
          :param int numY: number of divisions on y graph axis
          :param str interp: interpolation type when transferring to a grid. Valid values 'linear', 'nearest' or 'cubic'
          :param int fieldComponent: component of the field
          :param bool vertex: if vertices shoud be plot as points
          :param str colorBar: color bar details. Valid values '' for no colorbar, 'vertical' or 'horizontal'  
          :param str colorBarLegend: Legend for color bar. If '', current field name and units are printed. None prints nothing.
          :param tuple barRange: min and max bar range. If barRange=('NaN','NaN'), it is adjusted automatically
          :param str barFormatNum: format of color bar numbers
          :param str title: title
          :param str xlabel: x axis label
          :param str ylabel: y axis label
          :param str fileName: if nonempty, a filename is written to the disk, usually png, pdf, ps, eps and svg are supported
          :param bool show: if the plot should be showed
          :param tuple figsize: size of canvas in inches. Affects only showing a figure. Image to a file adjust one side automatically.
          :param obj matPlotFig: False means plot window remains in separate thread, True waits until a plot window becomes closed
          
          :return: handle to matPlotFig
          :rtype: matPlotFig
          """
          
          try:
              import numpy as np
              import math
              from scipy.interpolate import griddata
              import matplotlib
              matplotlib.use('TkAgg')  # Qt4Agg gives an empty, black window
              import matplotlib.pyplot as plt
          except ImportError as e:
              log.error('Skipping field2Image2D due to missing modules: %s' % e)
              return None
              # raise
          
          if self.fieldType != FieldType.FT_vertexBased:
              raise APIError.APIError('Only FieldType.FT_vertexBased is now supported')
          
          mesh = self.getMesh()
          numVertices = mesh.getNumberOfVertices()
  
          indX = 0
          indY = 0
          elev = 0
          if plane == 'xy':
              indX = 0
              indY = 1
              elev = 2
          elif plane == 'xz':
              indX = 0
              indY = 2
              elev = 1
          elif plane == 'yz':
              indX = 1
              indY = 2
              elev = 0
          
          # find eligible vertex points and values
          vertexPoints = []
          vertexValue = []
          for i in range(0, numVertices):
              coords = mesh.getVertex(i).getCoordinates()
              # print(coords)
              value = self.giveValue(i)[fieldComponent]
              
              if elevation[1] > coords[elev] > elevation[0]:
                  vertexPoints.append((coords[indX], coords[indY]))
                  vertexValue.append(value)
          
          if len(vertexPoints) == 0:
              log.info('No valid vertex points found, putting zeros on domain 1 x 1')
              for i in range(5):
                  vertexPoints.append((i % 2, i/4.))
                  vertexValue.append(0)
  
          # for i in range (0, len(vertexPoints)):
          #     print (vertexPoints[i], vertexValue[i])
  
          vertexPointsArr = np.array(vertexPoints)
          vertexValueArr = np.array(vertexValue)
          
          xMin = vertexPointsArr[:, 0].min()
          xMax = vertexPointsArr[:, 0].max()
          yMin = vertexPointsArr[:, 1].min()
          yMax = vertexPointsArr[:, 1].max()
          
          # print(xMin, xMax, yMin, yMax)
          
          grid_x, grid_y = np.mgrid[xMin:xMax:complex(0, numX), yMin:yMax:complex(0, numY)]
          grid_z1 = griddata(vertexPointsArr, vertexValueArr, (grid_x, grid_y), interp)
          
          # print (grid_z1.T)
          
          plt.ion()  # ineractive mode
          
          if matPlotFig is None:
              matPlotFig = plt.figure(figsize=figsize)
              # plt.xlim(xMin, xMax)
              # plt.ylim(yMin, yMax)
          
          plt.clf()
          plt.axis((xMin, xMax, yMin, yMax))
          image = plt.imshow(grid_z1.T, extent=(xMin, xMax, yMin, yMax), origin='lower', aspect='equal')
          # plt.margins(tight=True)
          # plt.tight_layout()
          # plt.margins(x=-0.3, y=-0.3)
  
          if colorBar:
              cbar = plt.colorbar(orientation=colorBar, format=barFormatNum)
              if colorBarLegend is not None:
                  if colorBarLegend == '':
                      colorBarLegend = self.getFieldIDName() + '_' + str(fieldComponent)
                      if self.unit is not None:
                          colorBarLegend = colorBarLegend + ' (' + self.unit.name() + ')'
                  cbar.set_label(colorBarLegend, rotation=0 if colorBar == 'horizontal' else 90)
          if title:
              plt.title(title)
          if xlabel:
              plt.xlabel(xlabel)
          if ylabel:
              plt.ylabel(ylabel)
          if vertex == 1:
              plt.scatter(vertexPointsArr[:, 0], vertexPointsArr[:, 1], marker='o', c='b', s=5, zorder=10)
  
          # plt.axis('equal')
          # plt.gca().set_aspect('equal', adjustable='box-forced')
          
          if isinstance(barRange[0], float) or isinstance(barRange[0], int):
              image.set_clim(vmin=barRange[0], vmax=barRange[1])
  
          if fileName:
              plt.savefig(fileName, bbox_inches='tight')
          if show:
              matPlotFig.canvas.draw()
              # plt.ioff()
              # plt.show(block=True)
          return matPlotFig
    
      def field2Image2DBlock(self):
          """
          Block an open window from matPlotLib. Waits until closed.
          """
          import matplotlib.pyplot as plt
          plt.ioff()
          plt.show(block=True)
  
      def toHdf5(self, fileName, group='component1/part1'):
          """
          Dump field to HDF5, in a simple format suitable for interoperability (TODO: document).
  
          :param str fileName: HDF5 file
          :param str group: HDF5 group the data will be saved under.
  
          The HDF hierarchy is like this::
  
              group
                |
                +--- mesh_01 {hash=25aa0aa04457}
                |      +--- [vertex_coords]
                |      +--- [cell_types]
                |      \--- [cell_vertices]
                +--- mesh_02 {hash=17809e2b86ea}
                |      +--- [vertex_coords]
                |      +--- [cell_types]
                |      \--- [cell_vertices]
                +--- ...
                +--- field_01
                |      +--- -> mesh_01
                |      \--- [vertex_values]
                +--- field_02
                |      +--- -> mesh_01
                |      \--- [vertex_values]
                +--- field_03
                |      +--- -> mesh_02
                |      \--- [cell_values]
                \--- ...
  
          where ``plain`` names are HDF (sub)groups, ``[bracketed]`` names are datasets, ``{name=value}`` are HDF attributes, ``->`` prefix indicated HDF5 hardlink (transparent to the user); numerical suffixes (``_01``, ...) are auto-allocated. Mesh objects are hardlinked using HDF5 hardlinks if an identical mesh is already stored in the group, based on hexdigest of its full data.
  
          .. note:: This method has not been tested yet. The format is subject to future changes.
          """
          import h5py
          hdf = h5py.File(fileName, 'a', libver='latest')
          if group not in hdf:
              gg = hdf.create_group(group)
          else:
              gg = hdf[group]
          # raise IOError('Path "%s" is already used in "%s".'%(path,fileName))
  
          def lowestUnused(trsf, predicate, start=1):
              """
              Find the lowest unused index, where *predicate* is used to test for existence, and *trsf* transforms
              integer (starting at *start* and incremented until unused value is found) to whatever predicate accepts
              as argument. Lowest transformed value is returned.
              """
              import itertools
              for i in itertools.count(start=start):
                  t = trsf(i)
                  if not predicate(t):
                      return t
          # save mesh (not saved if there already)
          newgrp = lowestUnused(trsf=lambda i: 'mesh_%02d' % i, predicate=lambda t: t in gg)
          mh5 = self.getMesh().asHdf5Object(parentgroup=gg, newgroup=newgrp)
  
          if self.value:
              fieldGrp = hdf.create_group(lowestUnused(trsf=lambda i, group=group: group+'/field_%02d' % i, predicate=lambda t: t in hdf))
              fieldGrp['mesh'] = mh5
              fieldGrp.attrs['fieldID'] = self.fieldID
              fieldGrp.attrs['valueType'] = self.valueType
              # string/bytes may not contain NULL when stored as string in HDF5
              # see http://docs.h5py.org/en/2.3/strings.html
              # that's why we cast to opaque type "void" and uncast using tostring before unpickling
              fieldGrp.attrs['units'] = numpy.void(pickle.dumps(self.unit))
              fieldGrp.attrs['time'] = numpy.void(pickle.dumps(self.time))
              # fieldGrp.attrs['time']=self.time.getValue()
              if self.fieldType == FieldType.FT_vertexBased:
                  val = numpy.empty(shape=(self.getMesh().getNumberOfVertices(), self.getRecordSize()), dtype=numpy.float)
                  for vert in range(self.getMesh().getNumberOfVertices()):
                      val[vert] = self.getVertexValue(vert).getValue()
                  fieldGrp['vertex_values'] = val
              elif self.fieldType == FieldType.FT_cellBased:
                  # raise NotImplementedError("Saving cell-based fields to HDF5 is not yet implemented.")
                  val = numpy.empty(shape=(self.getMesh().getNumberOfCells(), self.getRecordSize()), dtype=numpy.float)
                  for cell in range(self.getMesh().getNumberOfCells()):
                      val[cell] = self.getCellValue(cell)
                  fieldGrp['cell_values'] = val
              else:
                  raise RuntimeError("Unknown fieldType %d." % self.fieldType)
  
      @staticmethod
      def makeFromHdf5(fileName, group='component1/part1'):
          """
          Restore Fields from HDF5 file.
  
          :param str fileName: HDF5 file
          :param str group: HDF5 group the data will be read from (IOError is raised if the group does not exist).
          :return: list of new :obj:`Field` instances
          :rtype: [Field,Field,...]
  
  
          .. note:: This method has not been tested yet.
          """
          import h5py
          hdf = h5py.File(fileName, 'r', libver='latest')
          grp = hdf[group]
          # load mesh and field data from HDF5
          meshObjs = [obj for name, obj in grp.items() if name.startswith('mesh_')]
          fieldObjs = [obj for name, obj in grp.items() if name.startswith('field_')]
          # construct all meshes as mupif objects
          meshes = [Mesh.Mesh.makeFromHdf5Object(meshObj) for meshObj in meshObjs]
          # construct all fields as mupif objects
          ret = []
          for f in fieldObjs:
              if 'vertex_values' in f:
                  fieldType, values = FieldType.FT_vertexBased, f['vertex_values']
              elif 'cell_values' in f:
                  fieldType, values = FieldType.FT_cellBased, f['cell_values']
              else:
                  ValueError("HDF5/mupif format error: unable to determine field type.")
              fieldID, valueType, units, time = FieldID(f.attrs['fieldID']), f.attrs['valueType'], f.attrs['units'].tostring(), f.attrs['time'].tostring()
              if units == '':
                  units = None  # special case, handled at saving time
              else:
                  units = pickle.loads(units)
              if time == '':
                  time = None  # special case, handled at saving time
              else:
                  time = pickle.loads(time)
             
              meshIndex = meshObjs.index(f['mesh'])  # find which mesh object this field refers to
              ret.append(Field(mesh=meshes[meshIndex], fieldID=fieldID, units=units, time=time, valueType=valueType, values=values, fieldType=fieldType))
          return ret
  
      def toVTK2(self, fileName, format='ascii'):
          """
          Save the instance as Unstructured Grid in VTK2 format (``.vtk``).
  
          :param str fileName: where to save
          :param str format: one of ``ascii`` or ``binary``
          """
          self.field2VTKData().tofile(filename=fileName, format=format)
  
      @staticmethod
      def makeFromVTK2(fileName, unit, time=0, skip=['coolwarm']):
          """
          Return fields stored in *fileName* in the VTK2 (``.vtk``) format.
  
          :param str fileName: filename to load from
          :param PhysicalUnit unit: physical unit of filed values
          :param float time: time value for created fields (time is not saved in VTK2, thus cannot be recovered)
          :param [string,] skip: file names to be skipped when reading the input file; the default value skips the default coolwarm colormap.
          
          :returns: one field from VTK
          :rtype: Field 
          
          """
          import pyvtk
          from .dataID import FieldID
          if not fileName.endswith('.vtk'):
              log.warning('Field.makeFromVTK2: fileName should end with .vtk, you may get in trouble (proceeding).')
          ret = []
          try:
              data = pyvtk.VtkData(fileName)  # this is where reading the file happens (inside pyvtk)
          except NotImplementedError:
              log.info('pyvtk fails to open (binary?) file "%s", trying through vtk.vtkGenericDataReader.' % fileName)
              return Field.makeFromVTK3(fileName, time=time, units=unit, forceVersion2=True)
          ugr = data.structure
          if not isinstance(ugr, pyvtk.UnstructuredGrid):
              raise NotImplementedError(
                  "grid type %s is not handled by mupif (only UnstructuredGrid is)." % ugr.__class__.__name__)
          mesh = Mesh.UnstructuredMesh.makeFromPyvtkUnstructuredGrid(ugr)
          # get cell and point data
          pd, cd = data.point_data.data, data.cell_data.data
          for dd, fieldType in (pd, FieldType.FT_vertexBased), (cd, FieldType.FT_cellBased):
              for d in dd:
                  # will raise KeyError if fieldID with that name is not defined
                  if d.name in skip:
                      continue
                  fid = FieldID[d.name]
                  # determine the number of components using the expected number of values from the mesh
                  expectedNumVal = (mesh.getNumberOfVertices() if fieldType == FieldType.FT_vertexBased else mesh.getNumberOfCells())
                  nc = len(d.scalars)//expectedNumVal
                  valueType = ValueType.fromNumberOfComponents(nc)
                  values = [d.scalars[i*nc:i*nc+nc] for i in range(len(d.scalars))]
                  ret.append(Field(
                      mesh=mesh,
                      fieldID=fid,
                      units=unit,  # not stored at all
                      time=time,  # not stored either, set by caller
                      valueType=valueType,
                      values=values,
                      fieldType=fieldType
                  ))
          return ret
  
      def toVTK3(self, fileName, **kw):
          """
          Save the instance as Unstructured Grid in VTK3 format (``.vtu``). This is a simple proxy for calling :obj:`manyToVTK3` with the instance as the only field to be saved. If multiple fields with identical mesh are to be saved in VTK3, use :obj:`manyToVTK3` directly.
  
          :param fileName: output file name
          :param ``**kw``: passed to :obj:`manyToVTK3`
          """
          return self.manyToVTK3([self], fileName, **kw)
  
      @staticmethod
      def manyToVTK3(fields, fileName, ascii=False, compress=True):
          """
          Save all fields passed as argument into VTK3 Unstructured Grid file (``*.vtu``).
  
          All *fields* must be defined on the same mesh object; exception will be raised if this is not the case.
  
          :param list of Field fields:
          :param fileName: output file name
          :param bool ascii: write numbers are ASCII in the XML-based VTU file (rather than base64-encoded binary in XML)
          :param bool compress: apply compression to the data
          """
          import vtk
          if not fields:
              raise ValueError('At least one field must be passed.')
          # check if all fields are defined on the same mesh
          if len(set([f.mesh for f in fields])) != 1:
              raise RuntimeError(
                  'Not all fields are sharing the same Mesh object (and could not be saved to a single .vtu file')
          # convert mesh to VTK UnstructuredGrid
          mesh = fields[0].getMesh()
          vtkgrid = mesh.asVtkUnstructuredGrid()
          # add fields as arrays
          for f in fields:
              arr = vtk.vtkDoubleArray()
              arr.SetNumberOfComponents(f.getRecordSize())
              arr.SetName(f.getFieldIDName())
              assert f.getFieldType() in (FieldType.FT_vertexBased, FieldType.FT_cellBased)  # other future types not handled
              if f.getFieldType() == FieldType.FT_vertexBased:
                  nn = mesh.getNumberOfVertices()
              else:
                  nn = mesh.getNumberOfCells()
              arr.SetNumberOfValues(nn)
              for i in range(nn):
                  arr.SetTuple(i, f.giveValue(i))
              if f.getFieldType() == FieldType.FT_vertexBased:
                  vtkgrid.GetPointData().AddArray(arr)
              else:
                  vtkgrid.GetCellData().AddArray(arr)
          # write the unstructured grid to file
          writer = vtk.vtkXMLUnstructuredGridWriter()
          if compress:
              writer.SetCompressor(vtk.vtkZLibDataCompressor())
          if ascii:
              writer.SetDataModeToAscii()
          writer.SetFileName(fileName)
          # change between VTK5 and VTK6
          if vtk.vtkVersion().GetVTKMajorVersion() == 6:
              writer.SetInputData(vtkgrid)
          else:
              writer.SetInputData(vtkgrid)
          writer.Write()
          # finito
  
      @staticmethod
      def makeFromVTK3(fileName, units, time=0, forceVersion2=False):
          """
          Create fields from a VTK unstructured grid file (``.vtu``, format version 3, or ``.vtp`` with *forceVersion2*); the mesh is shared between fields.
  
          ``vtk.vtkXMLGenericDataObjectReader`` is used to open the file (unless *forceVersion2* is set), but it is checked that contained dataset is a ``vtk.vtkUnstructuredGrid`` and an error is raised if not.
  
          .. note:: Units are not supported when loading from VTK, all fields will have ``None`` unit assigned.
  
          :param str fileName: VTK (``*.vtu``) file
          :param PhysicalUnit units: units of read values
          :param float time: time value for created fields (time is not saved in VTK3, thus cannot be recovered)
          :param bool forceVersion2: if ``True``, ``vtk.vtkGenericDataObjectReader`` (for VTK version 2) will be used to open the file, isntead of ``vtk.vtkXMLGenericDataObjectReader``; this also supposes *fileName* ends with ``.vtk`` (not checked, but may cause an error).
          :return: list of new :obj:`Field` instances
          :rtype: [Field,Field,...]
          """
          import vtk
          from .dataID import FieldID
          # rr=vtk.vtkXMLUnstructuredGridReader()
          if forceVersion2 or fileName.endswith('.vtk'):
              rr = vtk.vtkGenericDataObjectReader()
          else:
              rr = vtk.vtkXMLGenericDataObjectReader()
          rr.SetFileName(fileName)
          rr.Update()
          ugrid = rr.GetOutput()
          if not isinstance(ugrid, vtk.vtkUnstructuredGrid):
              raise RuntimeError("vtkDataObject read from '%s' must be a vtkUnstructuredGrid (not a %s)" % (
                  fileName, ugrid.__class__.__name__))
          # import sys
          # sys.stderr.write(str((ugrid,ugrid.__class__,vtk.vtkUnstructuredGrid)))
          # make mesh -- implemented separately
          mesh = Mesh.UnstructuredMesh.makeFromVtkUnstructuredGrid(ugrid)
          # fields which will be returned
          ret = []
          # get cell and point data
          cd, pd = ugrid.GetCellData(), ugrid.GetPointData()
          for data, fieldType in (pd, FieldType.FT_vertexBased), (cd, FieldType.FT_cellBased):
              for idata in range(data.GetNumberOfArrays()):
                  aname, arr = pd.GetArrayName(idata), pd.GetArray(idata)
                  nt = arr.GetNumberOfTuples()
                  if nt == 0:
                      raise RuntimeError("Zero values in field '%s', unable to determine value type." % aname)
                  t0 = arr.GetTuple(0)
                  valueType = ValueType.fromNumberOfComponents(len(arr.GetTuple(0)))
                  # this will raise KeyError if fieldID with that name not defined
                  fid = FieldID[aname]
                  # get actual values as tuples
                  values = [arr.GetTuple(t) for t in range(nt)]
                  ret.append(Field(
                      mesh=mesh,
                      fieldID=fid,
                      units=units,  # not stored at all
                      time=time,  # not stored either, set by caller
                      valueType=valueType,
                      values=values,
                      fieldType=fieldType
                  ))
          return ret
  
      def _sum(self, other, sign1, sign2):
          """
          Should return a new instance. As deep copy is expensive,
          this operation should be avoided. Better to modify the field values.
          """
          raise TypeError('Not supported')
  
      def inUnitsOf(self, *units):
          """
          Should return a new instance. As deep copy is expensive,
          this operation should be avoided. Better to use convertToUnits method
          performing in place conversion.
          """
          raise TypeError('Not supported')
  
  #    def __deepcopy__(self, memo):
  #        """ Deepcopy operatin modified not to include attributes starting with underscore.
  #            These are supposed to be the ones valid only to s specific copy of the receiver.
  #            An example of these attributes are _PyroURI (injected by Application),
  #            where _PyroURI contains the URI of specific object, the copy should receive
  #            its own URI
  #        """
  #        cls = self.__class__
  #        dpcpy = cls.__new__(cls)
  #
  #        memo[id(self)] = dpcpy
  #        for attr in dir(self):
  #            if not attr.startswith('_'):
  #                value = getattr(self, attr)
  #                setattr(dpcpy, attr, copy.deepcopy(value, memo))
  #        return dpcpy
  