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

import sys
import vtk
import os

label = 'remesh_temp_'

def vtk2stl(infile, outfile):
    vtkReader = vtk.vtkUnstructuredGridReader()
    vtkReader.SetFileName(infile)
    vtkReader.Update()

    surfFilter = vtk.vtkDataSetSurfaceFilter()
    surfFilter.SetInputConnection(vtkReader.GetOutputPort())
    surfFilter.Update()

    stlWriter = vtk.vtkSTLWriter()
    stlWriter.SetFileName(outfile)
    stlWriter.SetInputConnection(surfFilter.GetOutputPort())
    stlWriter.Write()

geo_file = """
Mesh.RemeshAlgorithm=1;
Mesh.CharacteristicLengthFactor=$LEN;
Mesh.Algorithm3D = 1;

Merge "$FILE";
CreateTopology;

Compound Surface(100)={1};
//Physical Surface(101)={100};
Surface Loop(102)={100};

Volume(200)={102};
Physical Volume(201)={200};
"""

def remesh_create_geo(stlfile, geofile, length=0.7):

    f = open(geofile, 'wt')
    f.write(geo_file.replace('$LEN', str(length)).replace('$FILE', stlfile))
    f.close()

def main(argv):

    if argv[1][-4:] == '.vtk':
        stl_file = label + argv[1] + '.stl'
        vtk2stl(argv[1], stl_file)
    elif argv[1][-4:] == '.stl':
        stl_file = argv[1]
    else:
        raise(ValueError)

    geo_file = label + argv[1] + '.geo'
    remesh_create_geo(stl_file, geo_file, length=0.07)
    os.system("gmsh -3 %s -o out.vtk" % (geo_file))
    # os.system("rm remesh_temp_*")

if __name__ == '__main__':
    main(sys.argv)