# -*- coding: utf-8 -*-

##
# @author Pr Magoules HPC Research Group, CentraleSupelec, France
# @author G. Gbikpi-Benissan
# @date 2015-11-30, 2016-02-07, 2016-05-27, 2016-12-13
# @version 1.2
#
# @class Mesh
#

# -- Third-party modules
import numpy

##
# @brief A set of operations for handling mesh datasets.
#
# Currently:
# - File I/O (Legacy VTK Unstructured Grid)
# - Data fields (adding)
#
class Mesh ( object ) :

  # ----------------------------------------------------------------------------
  # -- CLASS ATTRIBUTES
  # ----------------------------------------------------------------------------

  # Error status
  SUCCESS = 0
  FAILURE = 1

  # VTK version
  VTK_VERSION = "4.2"

  # Class description
  CLASS_NAME = "Mesh"
  CLASS_AUTHOR = "G. G.-Benissan, MRG, CentraleSupelec, France"
  METHODS = """
  __init__ ( self )

  ReadFromFileVtk (
        self,
        file_name )
  ReadHeaderFromFileVtk (
        self,
        file_name )
  ReadUnstructGridFromFileVtkAscii (
        self,
        file_name )

  ReadFieldFromFileVtk (
        self,
        file_name,
        array_name )
  ReadFieldFromFileVtkAscii (
        self,
        file_name,
        array_name )

  WriteToFileVtk (
        self,
        file_name,
        file_type = "ASCII",
        dataset_type = "UNSTRUCTURED_GRID",
        title = "Generated by " + CLASS_AUTHOR + ".")
  WriteUnstructGridToFileVtkAscii (
        self,
        file_name,
        title = "Generated by " + CLASS_AUTHOR + ".")

  AllocPointData (
        self,
        numb_array )
  AllocCellData (
        self,
        numb_array )
  AddPointData (
        self,
        array_name,
        numb_column,
        array_val )
  AddCellData (
        self,
        array_name,
        numb_column,
        array_val )
  """

  # Data fields
  DEFAULT_MAX_NUMB_FIELD = 8

  # ----------------------------------------------------------------------------
  # -- INITIALIZATION
  # ----------------------------------------------------------------------------

  def __init__ ( self ) :

    # -- Dataset

    # Spatial dimension
    self.space_dim = 3

    # Number of nodes
    self.numb_node = 0

    # Nodes coordinates
    self.node_coord = numpy.array([])

    # Number of elements
    self.numb_elem = 0

    # Nodes of each element (CSR storage)
    # Values
    self.elem2node = numpy.array([])
    # Indices
    self.p_elem2node = numpy.array([0])

    # Elements types
    self.elem_type = numpy.array([])

    # -- Data fields

    # Number of data arrays
    self.numb_pdata = 0
    self.numb_cdata = 0

    # Data names
    self.pdata_name = Mesh.DEFAULT_MAX_NUMB_FIELD * [""]
    self.cdata_name = Mesh.DEFAULT_MAX_NUMB_FIELD * [""]

    # Data dimension (number of components of a value) in each array
    self.pdata_dim = numpy.array(Mesh.DEFAULT_MAX_NUMB_FIELD * [0])
    self.cdata_dim = numpy.array(Mesh.DEFAULT_MAX_NUMB_FIELD * [0])

    # Reference to data arrays
    self.pdata_val = Mesh.DEFAULT_MAX_NUMB_FIELD * [numpy.array([])]
    self.cdata_val = Mesh.DEFAULT_MAX_NUMB_FIELD * [numpy.array([])]

    # -- Error handling

    # Last error code
    self.err_code = Mesh.SUCCESS

    # Last error message
    self.err_msg = ""

  # END def __init__ ( self ) :
  # ----------------------------------------------------------------------------

  # ----------------------------------------------------------------------------
  # -- FILE I/O
  # ----------------------------------------------------------------------------

  ##
  # @brief Reads mesh dataset from a VTK file.
  # @param file_name = full name of the file (str)
  #
  def ReadFromFileVtk (
        self,
        file_name ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME + ".ReadFromFileVtk]"

    # -- read file header
    (file_type, dataset_type) = self.ReadHeaderFromFileVtk(file_name)
    if (self.err_code == Mesh.FAILURE) :
      self.err_msg = err_header + "\n" + self.err_msg
      return

    # -- read dataset
    if (dataset_type == "UNSTRUCTURED_GRID") :
      if (file_type == "ASCII") :

        self.ReadUnstructGridFromFileVtkAscii(file_name)
        if (self.err_code == Mesh.FAILURE) :
          self.err_msg = err_header + "\n" + self.err_msg

      else :
        self.err_code = Mesh.FAILURE
        self.err_msg = err_header + " Error: " + file_type
        self.err_msg += " file not supported for " + dataset_type + " dataset"

    else :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: " + dataset_type
      self.err_msg += " dataset not supported"


    return

  # END def ReadFromFileVtk (
#        self,
#        file_name ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Reads header part of a VTK file.
  # @param file_name = full name of the file (str)
  # @return file_type = numerical data format ("ASCII", "BINARY")
  # @return dataset_type = VTK dataset type ("UNSTRUCTURED_GRID", ...)
  # @return (file_type, dataset_type)
  #
  def ReadHeaderFromFileVtk (
        self,
        file_name ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME + ".ReadHeaderFromFileVtk]"

    # output
    file_type = ""
    dataset_type = ""

    # -- open file
    try :
      p_file = open(file_name, "r")
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot open " + file_name
      return (file_type, dataset_type)

    # -- skip version (# vtk DataFile Version x.y)
    buf = p_file.read(1)
    while (buf not in {'', '#'}) :
      buf = p_file.read(1)
    p_file.readline()

    # -- skip title
    p_file.readline()

    # -- read file type (ASCII|BINARY)

    # read first letter (A|B)
    buf = p_file.read(1)
    buf = buf.upper()
    while (buf not in {'', 'A', 'B'}) :
      buf = p_file.read(1)
      buf = buf.upper()
    file_type = buf

    # read remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      file_type = file_type + buf.upper()
      buf = p_file.read(1)

    # -- read dataset type (DATASET type)
    # type = STRUCTURED_POINTS|STRUCTURED_GRID|UNSTRUCTURED_GRID|POLYDATA
    #        |RECTILINEAR_GRID|FIELD

    # skip first letter (D)
    buf = p_file.read(1)
    buf = buf.upper()
    while (buf not in {'', 'D'}) :
      buf = p_file.read(1)
      buf = buf.upper()

    # skip remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # read type first letter (S|U|P|R|F)
    buf = p_file.read(1)
    buf = buf.upper()
    while (buf not in {'', 'S', 'U', 'P', 'R', 'F'}) :
      buf = p_file.read(1)
      buf = buf.upper()
    dataset_type = buf

    # read type remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      dataset_type = dataset_type + buf.upper()
      buf = p_file.read(1)

    # -- close file
    p_file.close()

    return (file_type, dataset_type)

  # END def ReadHeaderFromFileVtk (
#        self,
#        file_name ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Reads mesh dataset from a VTK "unstructured grid" ASCII file.
  # @remarks See ReadFromFileVtk
  #
  def ReadUnstructGridFromFileVtkAscii (
        self,
        file_name ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME
    err_header += ".ReadUnstructGridFromFileVtkAscii]"

    # output
    self.space_dim = 3
    self.numb_node = 0
    self.node_coord = numpy.array([])
    self.numb_elem = 0
    self.elem2node = numpy.array([])
    self.p_elem2node = numpy.array([0])
    self.elem_type = numpy.array([])

    # -- open file
    try :
      p_file = open(file_name, "r")
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot open " + file_name
      return

    # -- skip version (# vtk DataFile Version x.y)
    buf = p_file.read(1)
    while (buf not in {'', '#'}) :
      buf = p_file.read(1)
    p_file.readline()

    # -- skip title
    p_file.readline()

    # -- skip POINTS keyword (POINTS numb_node data_type)

    # skip first letter (P)
    buf = p_file.read(1)
    buf = buf.upper()
    while (buf not in {'', 'P'}) :
      buf = p_file.read(1)
      buf = buf.upper()

    # skip remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # -- read number of nodes

    # read first numeral symbol
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)
    self.numb_node = buf

    # read remaining numeral symbols
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      self.numb_node = self.numb_node + buf
      buf = p_file.read(1)

    # convert to number
    self.numb_node = int(self.numb_node)

    # -- skip data type

    # skip first letter
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)

    # skip remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # -- read nodes coordinates
    self.node_coord = []
    # for each node
    for i in range(0, self.numb_node) :
      self.node_coord.append([])
      # for each coordinate
      for j in range(0, self.space_dim) :

        # read first numeral symbol
        buf = p_file.read(1)
        while (buf in {' ', '\n', '\r'}) :
          buf = p_file.read(1)
        self.node_coord[i].append(buf)

        # read remaining numeral symbols
        buf = p_file.read(1)
        while (buf not in {'', ' ', '\n', '\r'}) :
          self.node_coord[i][j] = self.node_coord[i][j] + buf
          buf = p_file.read(1)

        # convert to number
        self.node_coord[i][j] = float(self.node_coord[i][j])

    # -- convert to numpy.ndarray
    self.node_coord = numpy.array(self.node_coord)

    # -- skip CELLS keyword (CELLS numb_elem size)

    # skip first letter (C)
    buf = p_file.read(1)
    buf = buf.upper()
    while (buf not in {'', 'C'}) :
      buf = p_file.read(1)
      buf = buf.upper()

    # skip remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # -- read number of elements

    # read first numeral symbol
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)
    self.numb_elem = buf

    # read remaining numeral symbols
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      self.numb_elem = self.numb_elem + buf
      buf = p_file.read(1)

    # convert to number
    self.numb_elem = int(self.numb_elem)

    # -- skip size

    # skip first numeral symbol
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)

    # skip remaining numeral symbols
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # -- read elements nodes
    self.elem2node = []
    self.p_elem2node = [0]
    # for each element
    for i in range(0, self.numb_elem) :

      # -- read number of nodes

      # read first numeral symbol
      buf = p_file.read(1)
      while (buf in {' ', '\n', '\r'}) :
        buf = p_file.read(1)
      elem_numb_node = buf

      # read remaining numeral symbols
      buf = p_file.read(1)
      while (buf not in {'', ' ', '\n', '\r'}) :
        elem_numb_node = elem_numb_node + buf
        buf = p_file.read(1)

      # convert to number
      elem_numb_node = int(elem_numb_node)

      # -- read nodes
      for j in range(0, elem_numb_node) :

        # read first numeral symbol
        buf = p_file.read(1)
        while (buf in {' ', '\n', '\r'}) :
          buf = p_file.read(1)
        self.elem2node.append(buf)

        # read remaining numeral symbols
        buf = p_file.read(1)
        while (buf not in {'', ' ', '\n', '\r'}) :
          self.elem2node[self.p_elem2node[i]+j] = self.elem2node[
                                                    self.p_elem2node[i]+j] + buf
          buf = p_file.read(1)

        # convert to number
        self.elem2node[self.p_elem2node[i]+j] = int(
                                          self.elem2node[self.p_elem2node[i]+j])

      # -- add index of next element
      self.p_elem2node.append(self.p_elem2node[i] + elem_numb_node)

    # -- convert to numpy.ndarray
    self.elem2node = numpy.array(self.elem2node)
    self.p_elem2node = numpy.array(self.p_elem2node)

    # -- skip CELL_TYPES keyword (CELL_TYPES numb_elem)

    # skip first letter (C)
    buf = p_file.read(1)
    buf = buf.upper()
    while (buf not in {'', 'C'}) :
      buf = p_file.read(1)
      buf = buf.upper()

    # skip remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # -- skip number of elements

    # skip first numeral symbol
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)

    # skip remaining numeral symbols
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n'}) :
      buf = p_file.read(1)

    # -- read elements types
    self.elem_type = []
    # for each element
    for i in range(0, self.numb_elem) :

      # -- read type

      # read first numeral symbol
      buf = p_file.read(1)
      while (buf in {' ', '\n', '\r'}) :
        buf = p_file.read(1)
      self.elem_type.append(buf)

      # read remaining numeral symbols
      buf = p_file.read(1)
      while (buf not in {'', ' ', '\n', '\r'}) :
        self.elem_type[i] = self.elem_type[i] + buf
        buf = p_file.read(1)

      # convert to number
      self.elem_type[i] = int(self.elem_type[i])

    # -- convert to numpy.ndarray
    self.elem_type = numpy.array(self.elem_type)

    # -- close file
    p_file.close()

    return

  # END def ReadUnstructGridFromFileVtkAscii (
#        self,
#        file_name ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Reads a field array from a VTK file.
  # @param file_name = full name of the file (str)
  # @param array_name = name of the array in the file (str)
  # @return array_data = tuples in the array (2D numpy.ndarray)
  #
  def ReadFieldFromFileVtk (
        self,
        file_name,
        array_name ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME + ".ReadFieldFromFileVtk]"

    # -- read file header
    (file_type, dataset_type) = self.ReadHeaderFromFileVtk(file_name)
    if (self.err_code == Mesh.FAILURE) :
      self.err_msg = err_header + "\n" + self.err_msg
      return

    # -- read field
    if (file_type == "ASCII") :

      array_data = self.ReadFieldFromFileVtkAscii(file_name, array_name)
      if (self.err_code == Mesh.FAILURE) :
        self.err_msg = err_header + "\n" + self.err_msg

    else :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: " + file_type
      self.err_msg += " file not supported for FIELD data"


    return array_data

  # END def ReadFieldFromFileVtk (
#        self,
#        file_name,
#        array_name ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Reads a field array from a VTK ASCII file.
  # @remarks See ReadFieldFromFileVtk
  #
  def ReadFieldFromFileVtkAscii (
        self,
        file_name,
        array_name ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME
    err_header += ".ReadFieldFromFileVtkAscii]"

    # output
    array_data = numpy.array([])

    # -- open file
    try :
      p_file = open(file_name, "r")
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot open " + file_name
      return array_data

    # -- skip version (# vtk DataFile Version x.y)
    buf = p_file.read(1)
    while (buf not in {'', '#'}) :
      buf = p_file.read(1)
    p_file.readline()

    # -- skip title
    p_file.readline()

    # -- read array name (arrayName numComponents numTuples dataType)

    buf2 = "#" + array_name
    while (buf2 not in {"", array_name}) :

      # read first letter of the array name
      buf1 = p_file.read(1)
      while (buf1 not in {'', array_name[0]}) :
        buf1 = p_file.read(1)
      buf2 = buf1

      # read remaining letters
      buf1 = p_file.read(1)
      while (buf1 not in {'', ' ', '\n', '\r'}) :
        buf2 += buf1
        buf1 = p_file.read(1)

    # -- read number of components (arrayName numComponents numTuples dataType)

    # read first numeral symbol
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)
    numb_component = buf

    # read remaining numeral symbols
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      numb_component += buf
      buf = p_file.read(1)

    # convert to number
    numb_component = int(numb_component)

    # -- read number of tuples (arrayName numComponents numTuples dataType)

    # read first numeral symbol
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)
    numb_tuple = buf

    # read remaining numeral symbols
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      numb_tuple += buf
      buf = p_file.read(1)

    # convert to number
    numb_tuple = int(numb_tuple)

    # -- read data type (arrayName numComponents numTuples dataType)

    # read first letter
    buf = p_file.read(1)
    while (buf in {' ', '\n', '\r'}) :
      buf = p_file.read(1)
    data_type = buf

    # read remaining letters
    buf = p_file.read(1)
    while (buf not in {'', ' ', '\n', '\r'}) :
      data_type += buf
      buf = p_file.read(1)

    # -- map VTK data type to Numpy data type
    data_type = data_type.lower()
    if (data_type in {"bit", "unsigned_short"}) :
      data_type = "short"
    elif (data_type in {"unsigned_char", "char"}) :
      data_type = "<U4"
    elif (data_type in {"unsigned_int", "vtkidtype"}) :
      data_type = "int"
    elif (data_type in {"unsigned_long"}) :
      data_type = "int64"

    # -- read data

    array_data = numpy.zeros((numb_tuple, numb_component), dtype=data_type)
    # for each tuple
    for i in range(0, numb_tuple) :
      # for each component
      for j in range(0, numb_component) :

        # read first character
        buf = p_file.read(1)
        while (buf in {' ', '\n', '\r'}) :
          buf = p_file.read(1)
        component = buf

        # read remaining characters
        buf = p_file.read(1)
        while (buf not in {'', ' ', '\n', '\r'}) :
          component += buf
          buf = p_file.read(1)

        # convert to data type
        array_data[i,j] = component # automatic conversion

    # -- close file
    p_file.close()

    return array_data

  # END def ReadFieldFromFileVtkAscii (
#        self,
#        file_name,
#        array_name ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Writes mesh dataset to a VTK file.
  # @param file_name = full name of the file (str)
  # @param file_type = numerical data format ("ASCII", "BINARY")
  #             [default: "ASCII"]
  # @param dataset_type = VTK dataset type ("UNSTRUCTURED_GRID", ...)
  #             [default: "UNSTRUCTURED_GRID"]
  # @param title = short description (255 char. max.)
  #             [default: "Generated by " + Mesh.CLASS_AUTHOR + "."]
  #
  def WriteToFileVtk (
        self,
        file_name,
        file_type = "ASCII",
        dataset_type = "UNSTRUCTURED_GRID",
        title = "Generated by " + CLASS_AUTHOR + "." ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME + ".WriteToFileVtk]"

    # -- write dataset
    file_type = file_type.upper()
    dataset_type = dataset_type.upper()
    if (dataset_type == "UNSTRUCTURED_GRID") :
      if (file_type == "ASCII") :
        self.WriteUnstructGridToFileVtkAscii(file_name, title)
        if (self.err_code == Mesh.FAILURE) :
          self.err_msg = err_header + "\n" + self.err_msg
      else :
        self.err_code = Mesh.FAILURE
        self.err_msg = err_header + "*** Error: " + file_type
        self.err_msg += " file not supported for " + dataset_type + " dataset"
    else :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + "*** Error: " + dataset_type
      self.err_msg += " dataset not supported"

    return

  # END def WriteToFileVtk (
#        self,
#        file_name,
#        file_type = "ASCII",
#        dataset_type = "UNSTRUCTURED_GRID",
#        title = "Generated by " + Mesh.CLASS_AUTHOR + ".") :
  # ----------------------------------------------------------------------------

  ##
  # @brief Writes mesh dataset to a VTK unstructured grid ASCII file.
  # @remarks See WriteToFileVtk
  #
  def WriteUnstructGridToFileVtkAscii (
        self,
        file_name,
        title = "Generated by " + CLASS_AUTHOR + "." ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME + ".WriteUnstructGridToFileVtkAscii]"

    # -- open file
    try :
      p_file = open(file_name, "w")
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot open " + file_name
      return

    # -- write version (# vtk DataFile Version x.y)
    p_file.write("# vtk DataFile Version " + Mesh.VTK_VERSION)

    # -- write title
    p_file.write("\n" + title + "\n")

    # -- write file type
    p_file.write("ASCII")

    # -- write dataset type
    p_file.write("\nDATASET UNSTRUCTURED_GRID")

    # -- write number of nodes
    p_file.write("\nPOINTS " + str(self.numb_node) + " double")

    # -- write nodes coordinates
    for i in range(0, self.numb_node) :
      p_file.write("\n" + str(self.node_coord[i][0]))
      p_file.write(" " + str(self.node_coord[i][1]))
      p_file.write(" " + str(self.node_coord[i][2]))

    # -- write number of elements
    p_file.write("\nCELLS " + str(self.numb_elem))

    # -- write elements list size
    p_file.write(" " + str(self.numb_elem + self.p_elem2node[self.numb_elem]))

    # -- write elements nodes
    # for each element
    for i in range(0, self.numb_elem) :

      # -- write number of nodes
      p_file.write("\n" + str(self.p_elem2node[i+1] - self.p_elem2node[i]))

      # -- write nodes
      for j in range(self.p_elem2node[i], self.p_elem2node[i+1]) :
        p_file.write(" " + str(self.elem2node[j]))

    # -- write number of elements
    p_file.write("\nCELL_TYPES " + str(self.numb_elem))

    # -- write elements types
    for i in range(self.numb_elem) :
      p_file.write("\n" + str(self.elem_type[i]))

    # -- write point data fields
    p_file.write("\nPOINT_DATA " + str(self.numb_node))
    p_file.write("\nFIELD AllPointData " + str(self.numb_pdata))
    for i in range(self.numb_pdata) :
      p_file.write("\n" + self.pdata_name[i] + " " + str(self.pdata_dim[i])
                   + " " + str(self.numb_node) + " double")
      for j in range(self.numb_node) :
        p_file.write("\n")
        for k in range(self.pdata_dim[i]) :
          p_file.write(str(self.pdata_val[i][j,k]) + " ")

    # -- write cell data fields
    p_file.write("\nCELL_DATA " + str(self.numb_elem))
    p_file.write("\nFIELD AllCellData " + str(self.numb_cdata))
    for i in range(self.numb_cdata) :
      p_file.write("\n" + self.cdata_name[i] + " " + str(self.cdata_dim[i])
                   + " " + str(self.numb_elem) + " double")
      for j in range(self.numb_elem) :
        p_file.write("\n")
        for k in range(self.cdata_dim[i]) :
          p_file.write(str(self.cdata_val[i][j,k]) + " ")

    # -- close file
    p_file.close()

    return

  # END def WriteUnstructGridToFileVtkAscii (
#        self,
#        file_name,
#        title = "Generated by " + CLASS_AUTHOR + ".") :
  # ----------------------------------------------------------------------------

  # ----------------------------------------------------------------------------
  # -- DATA FIELDS
  # ----------------------------------------------------------------------------

  ##
  # @brief Allocate memory for storing references to point data arrays.
  # @param numb_array = number of arrays
  #
  def AllocPointData (
        self,
        numb_array ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME
    err_header += ".AllocPointData]"

    # output
    self.numb_pdata = 0

    # -- allocate
    try :
      self.pdata_dim = numpy.array(numb_array * [0])
      self.pdata_name = numb_array * [""]
      self.pdata_val = numb_array * [numpy.array([])]
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot allocate memory for"
      self.err_msg += " referencing " + str(numb_array) + " arrays."
      return

    return

  # END def AllocPointData (
#        self,
#        numb_array ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Allocate memory for storing references to cell data arrays.
  # @param numb_array = number of arrays
  #
  def AllocCellData (
        self,
        numb_array ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME
    err_header += ".AllocCellData]"

    # output
    self.numb_cdata = 0

    # -- allocate
    try :
      self.cdata_dim = numpy.array(numb_array * [0])
      self.cdata_name = numb_array * [""]
      self.cdata_val = numb_array * [numpy.array([])]
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot allocate memory for"
      self.err_msg += " referencing " + str(numb_array) + " arrays."
      return

    return

  # END def AllocCellData (
#        self,
#        numb_array ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Add a reference to a point data array.
  # @param array_name = name of the array (for VTK file format)
  # @param numb_column = number of components of each value
  # @param array_val = reference to the array contents (2D numpy.ndarray)
  #
  def AddPointData (
        self,
        array_name,
        numb_column,
        array_val ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME
    err_header += ".AddPointData]"

    # -- add reference to array
    try :
      self.pdata_dim[self.numb_pdata] = numb_column
      self.pdata_name[self.numb_pdata] = array_name
      self.pdata_val[self.numb_pdata] = array_val
      self.numb_pdata += 1
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot add reference to "
      self.err_msg += str(array_name) + " array."
      return

    return

  # END def AddPointData (
#        self,
#        array_name,
#        numb_column,
#        array_val ) :
  # ----------------------------------------------------------------------------

  ##
  # @brief Add a reference to a cell data array.
  # @param array_name = name of the array (for VTK file format)
  # @param numb_column = number of components of each value
  # @param array_val = reference to the array contents (2D numpy.ndarray)
  #
  def AddCellData (
        self,
        array_name,
        numb_column,
        array_val ) :

    # -- init

    # error handling
    self.err_code = Mesh.SUCCESS
    self.err_msg = ""
    err_header = "*** [" + Mesh.CLASS_NAME
    err_header += ".AddCellData]"

    # -- add reference to array
    try :
      self.cdata_dim[self.numb_cdata] = numb_column
      self.cdata_name[self.numb_cdata] = array_name
      self.cdata_val[self.numb_cdata] = array_val
      self.numb_cdata += 1
    except :
      self.err_code = Mesh.FAILURE
      self.err_msg = err_header + " Error: cannot add reference to "
      self.err_msg += str(array_name) + " array."
      return

    return

  # END def AddCellData (
#        self,
#        array_name,
#        numb_column,
#        array_val ) :
  # ----------------------------------------------------------------------------

# END class Mesh ( object ) :
# ------------------------------------------------------------------------------