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
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
