You can use h5tovtk, which is part of the h5utils package to convert the
epsilon.h5 file into a VTK file which can be visualized in paraview or
mayavi.

If you want more periods, you can use "mpb-data" to generate them.
See here for an example:
http://ab-initio.mit.edu/wiki/index.php/MPB_Data_Analysis_Tutorial#Visualizing_the_diamond_lattice_structure_and_bands

The example uses h5tov5d and vis5d, but you can easily use h5tovtk and
paraview instead.

One problem with h5tovtk is that it does not take into account the
direction of the lattice vectors.
If you want to visualize structures based on non-cartesian lattices, you
can use my h5tovts.py script, which creates a VTK structured grid as
output with axes based on the lattice vectors.
cf attachments.

I haven't used these scripts much for a while, so some things may be
broken. Let me know if you have any problems.
It also does not yet create multiple periods like mpb-data
unfortunately, and runs a bit slow. A C/C++ version would probably be
faster.

This could also interest you:
https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg04989.html

Regards,
Mike


On 22/08/15 14:29, vijaykumar gudelli wrote:
> Dear MPB-developers and users,
> 
> I'm very new to MPB. I have succeeded in reproducing the tutorials and
> some of the examples too.
> Now I'm facing problem in the visualization of .h5 files with some test
> case structure in 3D.  The structure file which I would like to generate
> is similar to the Figure 5, of chapter 3 (page no. 33) "Photonic
> crystals Molding the flow of light" (Second edition). The problem here
> I'm facing is, how to visualize the file in the 3D? Is there any
> visualization tools available to view in 3D of file epsilon.h5?
> I found BEAM as one of the solution for it, here I can able to see all
> the files in the epsilon.h5 individually, but I couldn't get much of its
> information about the 3D structure.
> Please direct me to get the solution for this problem.
> Any help in the regards is highly appreciated.
> 
> Thanking you,
> -- 
> -- _Best_ _Regards_
> __
> V. K. GUDELLI
> vkgude...@gmail.com <mailto:vkgude...@gmail.com>
> 
> 
> _______________________________________________
> mpb-discuss mailing list
> mpb-discuss@ab-initio.mit.edu
> http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/mpb-discuss
> 


#!/usr/bin/env python2
# -*- coding: utf-8 -*-

import os
import sys
import argparse

from h5utils import h5tovts

def h5tovts_argparse():

  parser = argparse.ArgumentParser(description = 'Convert an MPB-created HDF5 file to a .vts file.')
  parser.add_argument('-v','--verbose', action="count", dest="verbosity", default=0, help='verbosity level')

  parser.add_argument('h5file', action="store", help='input HDF5 file', nargs='+')
  #parser.add_argument('vtsfile', action="store", help='output .vts file', nargs='?')
  #parser.add_argument('-b','--basepath', action="store", default=None, help='basepath for output files')

  parser.add_argument('--size', nargs=3, type=float, default=None, help='lattice size (only used for .h5 files without lattice vectors, i.e. MEEP output, not for MPB output.)')
  parser.add_argument('-d', '--dataset', help='use dataset <name> in the input files (default: first dataset)', metavar='name')

  arguments = parser.parse_args()

  if arguments.verbosity > 0:
    print('---------')
    print(arguments)
    print('---------')
    
  if not len(sys.argv) > 1:
    parser.print_help()
  else:
    #if arguments.basepath is None:
      #arguments.basepath = os.path.splitext(arguments.h5file)[0]

    for h5file in arguments.h5file:
      basepath = os.path.splitext(h5file)[0]
      print('==> {} -> {}.vts and {}.vti'.format(h5file, basepath, basepath))
      h5tovts(h5file, basepath, arguments.size, arguments.dataset, arguments.verbosity)
  
  return

if __name__ == '__main__':
  h5tovts_argparse()
#!/usr/bin/env python2
# -*- coding: utf-8 -*-

'''
This module provides various utilities to work with the HDF5 format (and the VTK formats).

Automatic documentation cannot be generated at the moment, since we use python3 by default.

This script only runs with python2 because the python VTK module is not yet fully ported to python3.

.. note:: VTK numbering is in increasing X, then Y, then Z.

.. todo:: Use virtual void vtkSTLReader::MergingOff() 	[virtual] - Turn on/off merging of points/triangles.
.. todo:: Read in multiple separate STL files.
'''

from __future__ import division

import os
import sys
import vtk
import h5py
import time
import numpy
from numpy import array, zeros, sqrt, linspace
from numpy.linalg import norm
from vtk.util.numpy_support import numpy_to_vtk

class Lattice(object):
  '''
  The lattice class is normally used only for the geometry-lattice variable, and specifies the three lattice directions of the crystal and the lengths of the corresponding lattice vectors.

  lattice
      Properties: 

  basis1, basis2, basis3 [vector3]
      The three lattice directions of the crystal, specified in the cartesian basis. The lengths of these vectors are ignored--only their directions matter. The lengths are determined by the basis-size property, below. These vectors are then used as a basis for all other 3-vectors in the ctl file. They default to the x, y, and z directions, respectively. 

  basis-size [vector3]
      The components of basis-size are the lengths of the three basis vectors, respectively. They default to unit lengths. 

  size [vector3]
      The size of the lattice (i.e. the length of the lattice vectors Ri, in which the crystal is periodic) in units of the basis vectors. Thus, the actual lengths of the lattice vectors are given by the components of size multiplied by the components of basis-size. (Alternatively, you can think of size as the vector between opposite corners of the primitive cell, specified in the lattice basis.) Defaults to unit lengths. 

  resolution [number or vector3]
      Specifies the computational grid resolution, in pixels per lattice unit (a lattice unit is one basis vector in a given direction).
      If resolution is a vector3, then specifies a different resolution for each direction;
      otherwise the resolution is uniform.
      (The grid size is then the product of the lattice size and the resolution, rounded up to the next positive integer.)
      Defaults to 10.
    
  If any dimension has the special size no-size, then the dimensionality of the problem is reduced by one; strictly speaking, the dielectric function is taken to be uniform along that dimension. (In this case, the no-size dimension should generally be orthogonal to the other dimensions.) 
    
  .. todo:: This class still needs some work. But works to currently store necessary variables. Use with care.
  .. todo:: merge somehow with BFDTD meshing classes? (need proper definition for resolution, etc)
  '''
  def __init__(self):
    self.basis1 = array([1,0,0])
    self.basis2 = array([0,1,0])
    self.basis3 = array([0,0,1])
    self.basis_size = array([1,1,1])
    self.size = array([1,1,1])
    
    self.xmesh = [0,1]
    self.ymesh = [0,1]
    self.zmesh = [0,1]
    self.setResolution(10, 10, 10)
    
    return
  
  def __str__(self):
    ret = 'lattice:\n'
    ret += '  basis1 = {}\n'.format(self.basis1)
    ret += '  basis2 = {}\n'.format(self.basis2)
    ret += '  basis3 = {}\n'.format(self.basis3)
    ret += '  basis_size = {}\n'.format(self.basis_size)
    ret += '  size = {}'.format(self.size)
    return(ret)
  
  def getLatticeVectors(self):
    a1 = self.size[0]*self.basis_size[0]*self.basis1/norm(self.basis1)
    a2 = self.size[1]*self.basis_size[1]*self.basis2/norm(self.basis2)
    a3 = self.size[2]*self.basis_size[2]*self.basis3/norm(self.basis3)
    return (a1, a2, a3)
  
  def getBounds(self):
    (a1, a2, a3) = self.getLatticeVectors()

    Pmax = 0.5*a1 + 0.5*a2 + 0.5*a3
    Pmin = -Pmax

    xmin = Pmin[0]
    ymin = Pmin[1]
    zmin = Pmin[2]

    xmax = Pmax[0]
    ymax = Pmax[1]
    zmax = Pmax[2]
    
    return (xmin,xmax, ymin,ymax, zmin,zmax)
    
  def getXmeshDelta(self):
    return(numpy.diff(self.xmesh))
  def getYmeshDelta(self):
    return(numpy.diff(self.ymesh))
  def getZmeshDelta(self):
    return(numpy.diff(self.zmesh))
  def getMeshDelta(self):
    return(numpy.diff(self.xmesh),numpy.diff(self.ymesh),numpy.diff(self.zmesh))
    
  def getMinDeltas(self):
    dx = min(self.getXmeshDelta())
    dy = min(self.getYmeshDelta())
    dz = min(self.getZmeshDelta())
    return (dx,dy,dz)

  def getResolution(self):
    return [len(self.xmesh), len(self.ymesh), len(self.zmesh)]

  def getSpacing(self):
    (xmin,xmax, ymin,ymax, zmin,zmax) = self.getBounds()
    (Nx, Ny, Nz) = self.getResolution()
    dx = (xmax-xmin)/(Nx-1)
    dy = (ymax-ymin)/(Ny-1)
    dz = (zmax-zmin)/(Nz-1)
    return (dx,dy,dz)
        
  def setResolution(self, Nx, Ny, Nz):
    '''
    Sets up homogeneous X,Y,Z grids with Nx,Ny,Nz points.
    '''
    self.xmesh = linspace(-0.5, 0.5, Nx)
    self.ymesh = linspace(-0.5, 0.5, Ny)
    self.zmesh = linspace(-0.5, 0.5, Nz)
    return

  def getMesh(self):
    return (self.xmesh, self.ymesh, self.zmesh)

  def setXmesh(self, xmesh):
    self.xmesh = xmesh
    return

  def setYmesh(self, ymesh):
    self.ymesh = ymesh
    return

  def setZmesh(self, zmesh):
    self.zmesh = zmesh
    return

  def setSize(self, size):
    self.size = size

  def setBasisSize(self, basis_size):
    self.basis_size = basis_size

class FCClattice(Lattice):
  ''' Create a FCC lattice. '''
  def __init__(self):
    ''' Constructor '''
    super(FCClattice, self).__init__()
    # set up lattice
    self.basis1 = array([0, 1, 1])
    self.basis2 = array([1, 0, 1])
    self.basis3 = array([1, 1, 0])
    self.basis_size = array([sqrt(0.5), sqrt(0.5), sqrt(0.5)])

class BCClattice(Lattice):
  ''' Create a BCC lattice. '''
  def __init__(self):
    ''' Constructor '''
    super(BCClattice, self).__init__()
    # set up lattice
    self.basis1 = array([-1,  1,  1])
    self.basis2 = array([ 1, -1,  1])
    self.basis3 = array([ 1,  1, -1])
    self.basis_size = array([sqrt(3)/2, sqrt(3)/2, sqrt(3)/2])
    
def h5tovts(h5file, basepath, total_lattice_size=None, dataset=None, verbosity=0):
  '''
  * total_lattice_size : total length of the lattice vectors, i.e. **size*basis-size** in MPB terms. Overrides any values obtained from reading the lattice vectors in the .h5 file.
  
  .. todo:: Emulate the -x/y/z options of mpb-data, i.e. create a periodic structure from a unit-cell.
  '''
  # read in .h5 file
  with h5py.File(h5file, "r") as f:
    print('Reading from ' + h5file)
    
    # get description
    description = None
    if 'description' in f.keys():
      #description              Dataset {SCALAR}
      description = f['description'][...].tostring().decode("ascii").strip('\0')
      print('description = {}'.format(description))

    # create list of datasets
    dataset_list = [k for k in f.keys() if k not in ['description','lattice vectors', 'xmesh', 'ymesh', 'zmesh'] ]
    print('Available datasets: {}'.format(dataset_list))

    # choose first dataset if not specified
    if dataset is None:
      dataset = dataset_list[0]

    # set up data, Nx, Ny, Nz
    # TODO: Might cause conflict with size read from x/y/z mesh values... Add warnings?
    print('Using dataset = {}'.format(dataset))
    data = f[dataset]
    (Nx, Ny, Nz) = data.shape
    
    # set up lattice
    mylattice = Lattice()
    
    if 'lattice vectors' in f.keys():
      #lattice\ vectors         Dataset {3, 3}
      lattice_vectors = f['lattice vectors'][...]
      mylattice.basis1 = lattice_vectors[0]
      mylattice.basis2 = lattice_vectors[1]
      mylattice.basis3 = lattice_vectors[2]
      mylattice.setSize( [norm(mylattice.basis1), norm(mylattice.basis2), norm(mylattice.basis3)] )

    if total_lattice_size:
      mylattice.setSize(total_lattice_size)
    
    # generate homogeneous mesh by default
    mylattice.setResolution(Nx, Ny, Nz)

    if 'xmesh' in f.keys():
      mylattice.setXmesh(f['xmesh'][...])
    if 'ymesh' in f.keys():
      mylattice.setYmesh(f['ymesh'][...])
    if 'zmesh' in f.keys():
      mylattice.setZmesh(f['zmesh'][...])
    
    (a1, a2, a3) = mylattice.getLatticeVectors()
    (xmesh, ymesh, zmesh) = mylattice.getMesh()
    print('a1 = {}'.format(a1))
    print('a2 = {}'.format(a2))
    print('a3 = {}'.format(a3))
    print('data.shape = {} x {} x {}'.format(Nx, Ny, Nz))
    print('mesh.shape = {} x {} x {}'.format(len(xmesh), len(ymesh), len(zmesh)))
    if Nx != len(xmesh) or Ny != len(ymesh) or Nz != len(zmesh):
      raise Exception('Inconsistent number of cells between the dataset and the x/y/z meshs.')
    
    # create the vtkPoints structure for the coordinates
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(Nx*Ny*Nz)

    # create the vtkFloatArray structures for the data
    dataset_dict = dict()
    for key in dataset_list:
      print('key = {}'.format(key))
      scalar = f[key][...].transpose().reshape(-1,1)
      vtk_data = numpy_to_vtk(scalar)
      #vtk_data = vtk.vtkFloatArray()
      vtk_data.SetName(key)
      #vtk_data.SetNumberOfTuples(Nx*Ny*Nz)
      dataset_dict[key] = (f[key], vtk_data, scalar)

    last_info_time = time.time()
    
    print('Starting loops')
    counter = 0
    # fill the vtkPoints and vtkFloatArray
    for k in range(Nz):
      for j in range(Ny):
        for i in range(Nx):
          offset = i + j*Nx + k*Nx*Ny

          # old system:
          #   coord = (i/(Nx-1) - 0.5)*a1 + (j/(Ny-1) - 0.5)*a2 + (k/(Nz-1) - 0.5)*a3
          # new system:
          coord = xmesh[i]*a1 + ymesh[j]*a2 + zmesh[k]*a3

          points.SetPoint(offset, coord)
          #for key in dataset_dict.keys():
            #dataset_dict[key][1].SetTuple1(offset, dataset_dict[key][0][i,j,k])
          
          #InsertTuples 	
          #virtual void vtkAbstractArray::InsertTuples 	( 	vtkIdList *  	dstIds,
          #vtkIdList *  	srcIds,
          #vtkAbstractArray *  	source 
          #) 		[pure virtual]

          if time.time() - last_info_time > 5:
            print('{} %'.format(100*offset/(Nx*Ny*Nz-1)))
            last_info_time = time.time()
  
          if verbosity>1:
            counter += 1
            progress_str = 'Progress: {}/{}'.format(counter, Nx*Ny*Nz)
            #print(progress_str, end='\r')
            print(progress_str)
            #subprocess.call(["printf", progress_str+'\r'])
    
    print('\nLoops done.')

    #for key in dataset_dict.keys():
      #h5_data = dataset_dict[key][0]
      #vtk_data = dataset_dict[key][1]
      #vtk_data.InsertTuples()
      #.SetTuple1(offset, [i,j,k])
    
    # create structured grid
    dataset_vts = vtk.vtkStructuredGrid()
    dataset_vts.SetDimensions(Nx, Ny, Nz)
    dataset_vts.SetPoints(points)
    
    # create vtkImageData
    dataset_vti = vtk.vtkImageData()
    dataset_vti.SetDimensions(Nx, Ny, Nz)
    (xmin,xmax, ymin,ymax, zmin,zmax) = mylattice.getBounds()
    dataset_vti.SetOrigin([xmin, ymin, zmin])
    dataset_vti.SetSpacing(mylattice.getSpacing())
    
    # add scalar data to the grids
    for key in dataset_dict.keys():
      dataset_vts.GetPointData().AddArray(dataset_dict[key][1])
      dataset_vti.GetPointData().AddArray(dataset_dict[key][1])
    
    dataset_vts.GetPointData().SetActiveScalars('data')
    dataset_vti.GetPointData().SetActiveScalars('data')

    # write out .vts file
    writer = vtk.vtkXMLStructuredGridWriter()
    writer.SetInputData(dataset_vts)
    writer.SetFileName(basepath + '.' + writer.GetDefaultFileExtension())
    writer.Write()

    # write out .vti file
    writer = vtk.vtkXMLImageDataWriter()
    writer.SetInputData(dataset_vti)
    writer.SetFileName(basepath + '.' + writer.GetDefaultFileExtension())
    writer.Write()

    return

def stltoh5(stlfile, basepath, epsilon_inside, epsilon_outside, lattice=Lattice(), verbosity=0):
  
  if not os.path.exists(stlfile):
    if sys.version_info.major == 2:
      raise IOError('No such file or directory: {}'.format(stlfile)) # py2
    else:
      raise FileNotFoundError('No such file or directory: {}'.format(stlfile)) # py3
  
  print('--> timer start')
  time_start = time.time()

  # read in .stl file
  reader = vtk.vtkSTLReader()
  reader.SetFileName(stlfile)
  reader.Update()
  polydata = reader.GetOutput()

  # write .vtp file
  writer = vtk.vtkXMLPolyDataWriter()
  writer.SetInputData(polydata)
  writer.SetFileName(basepath + '.' + writer.GetDefaultFileExtension())
  writer.Write()

  # set up implicit_function
  implicit_function = vtk.vtkImplicitPolyDataDistance()
  implicit_function.SetInput(polydata)

  print("--> Elapsed time: %.4f sec" % (time.time() - time_start))
  
  stl_to_vts_and_h5(implicit_function, basepath, lattice, epsilon_inside, epsilon_outside)
  print("--> Elapsed time: %.4f sec" % (time.time() - time_start))
  return

def stl_to_vts_and_h5(implicit_function, outfile_basename, lattice, epsilon_inside, epsilon_outside):
  (Nx, Ny, Nz) = lattice.getResolution()
  (a1, a2, a3) = lattice.getLatticeVectors()
  
  points = vtk.vtkPoints()
  points.SetNumberOfPoints(Nx*Ny*Nz)

  scalars_vtk = vtk.vtkFloatArray()
  scalars_vtk.SetNumberOfTuples(Nx*Ny*Nz)
  scalars_numpy = zeros([Nx,Ny,Nz])

  print('=== dims ===')
  print(scalars_numpy.shape)
  
  last_info_time = time.time()

  print('=== Loop start ===')
  for k in range(Nz):
    for j in range(Ny):
      for i in range(Nx):
        coord = (i/(Nx-1) - 0.5)*a1 + (j/(Ny-1) - 0.5)*a2 + (k/(Nz-1) - 0.5)*a3
        offset = i + j*Nx + k*Nx*Ny
        points.SetPoint(offset, coord)
        if implicit_function.FunctionValue(coord) <= 0:
          value = epsilon_inside
        else:
          value = epsilon_outside
        scalars_vtk.SetTuple1(offset, value)
        scalars_numpy[i, j, k] = value
        
        if time.time() - last_info_time > 5:
          print('{} %'.format(100*offset/(Nx*Ny*Nz-1)))
          last_info_time = time.time()
  
  print('=== Loop end ===')

  dataset = vtk.vtkStructuredGrid()
  dataset.SetDimensions(Nx, Ny, Nz)
  dataset.SetPoints(points)
  dataset.GetPointData().SetScalars(scalars_vtk)

  writer = vtk.vtkXMLStructuredGridWriter()
  writer.SetInputData(dataset)
  writer.SetFileName(outfile_basename + '.' + writer.GetDefaultFileExtension())
  writer.Write()
  
  h5file = outfile_basename + '.h5'
  with h5py.File(h5file, "w") as f:
    print('writing to ' + h5file)
    
    dset = f.create_dataset('/data', scalars_numpy.shape, dtype=numpy.float64)
    dset[...] = scalars_numpy

    dset = f.create_dataset("description", (), dtype="S29")
    dset[...] = 'dielectric function, epsilon'
    
    lattice_vectors = numpy.array([a1, a2, a3])
    print(lattice_vectors)

    dset = f.create_dataset('/lattice vectors', lattice_vectors.shape, dtype=numpy.float64)
    dset[...] = lattice_vectors
    
    # TODO: Add these fields:
    #epsilon.xx               Dataset {100, 100, 100}
    #epsilon.xy               Dataset {100, 100, 100}
    #epsilon.xz               Dataset {100, 100, 100}
    #epsilon.yy               Dataset {100, 100, 100}
    #epsilon.yz               Dataset {100, 100, 100}
    #epsilon.zz               Dataset {100, 100, 100}
    #epsilon_inverse.xx       Dataset {100, 100, 100}
    #epsilon_inverse.xy       Dataset {100, 100, 100}
    #epsilon_inverse.xz       Dataset {100, 100, 100}
    #epsilon_inverse.yy       Dataset {100, 100, 100}
    #epsilon_inverse.yz       Dataset {100, 100, 100}
    #epsilon_inverse.zz       Dataset {100, 100, 100}
  
  return

if __name__ == '__main__':
	pass

Attachment: signature.asc
Description: OpenPGP digital signature

_______________________________________________
mpb-discuss mailing list
mpb-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/mpb-discuss

Reply via email to