[cig-commits] commit: updated MADDs5 vtk script

[Mercurial]

changeset:   106:18d756141ac1
user:        Marc Spiegelman <mspieg at ldeo.columbia.edu>
date:        Mon Dec 07 00:03:09 2009 -0500
files:       MADDs-5/MADDs-5a/run/MADDs-5atovtk.py
description:
updated MADDs5 vtk script


diff -r 81e5071e5e51 -r 18d756141ac1 MADDs-5/MADDs-5a/run/MADDs-5atovtk.py
--- a/MADDs-5/MADDs-5a/run/MADDs-5atovtk.py	Sun Dec 06 21:54:16 2009 -0500
+++ b/MADDs-5/MADDs-5a/run/MADDs-5atovtk.py	Mon Dec 07 00:03:09 2009 -0500
@@ -8,8 +8,10 @@ from glob import glob
 from glob import glob
 from optparse import OptionParser
 
-# import vtk utility functions
-sys.path[0] += "/../../../python"
+#import vtk utility functions
+REPOHOME = subprocess.Popen(["hg","root"],stdout=subprocess.PIPE).communicate()[0]
+REPOHOME = REPOHOME.replace("\n","")
+sys.path.insert(0,REPOHOME+"/python")
 
 from vtkutils import *
 
@@ -32,10 +34,9 @@ else:
 
 # get arguments from directory name and input file
 
-# determine if using P1 or P2 elements
-ONAME = os.path.basename(ODIR)
-use_P2 = ONAME.find('P2') >= 0
-    
+# determine if using P1 or P2 elements by parsing the directory name 
+use_P2 = ODIR.find('P2') >= 0
+
 # get the mesh file
 iFile = open("input.options","r")
 s = iFile.readlines()
@@ -55,43 +56,52 @@ dim = mesh.topology().dim()
 
 
 #if using quadratic elements make a new finer mesh for P1
+P2 = FunctionSpace(mesh,"CG",2);
+P1 = FunctionSpace(mesh,"CG",1);
 if use_P2:
-    P2 = FunctionSpace(mesh,"CG",2);
     fine_mesh = Mesh(mesh)
     fine_mesh.refine()
-    P1 = FunctionSpace(fine_mesh, "CG",1)
-    V = P1  # this isn't correct but it will work
+    P1f = FunctionSpace(fine_mesh, "CG",1)
+    V = P1f  
 else:
-    P1 = FunctionSpace(mesh, "CG",1)
     V = P1
 
 # need a check for this
-M = V + V;
+M1 = V + V;
+M2 = P2 + P2;
 
 # get output file names (and compress if necessary)
 vector_files =  glob("u_vector*.xml")
+vector_files.sort()
 if len(vector_files) > 0:
     print 'compressing .xml files'
     subprocess.call(['gzip','-fv']+vector_files)
 
 vector_files =  glob("u_vector*.xml.gz")
+vector_files.sort()
 if len(vector_files) == 0:
     print 'Error: no compressed vector files'
     exit()
 
 pvd_file = "PressurePorosity.pvd"
-ufile = File(pvd_file)
+ufile = File(pvd_file,"compressed")
 for i in vector_files[:]:
-    u = Function(M,i) # this really shouldn't work but I'm hacking at the moment
-    ufile << u
-
+    if use_P2:
+        u = Function(M2,i) # this really shouldn't work but I'm hacking at the moment
+        u_out = Function(M1)
+        u_out.interpolate(u);
+        ufile << u_out
+    else:
+        u = Function(M1,i)
+        ufile << u
+        
 # get list of times
 t = []
 for i in vector_files[:]:
     s = i.partition('_t')[2]
     t.append(s.split('_')[0])
 
-# clean up pvd and vtk files to have proper times and names
+# clean up pvd and vtu files to have proper times and names
 pvd_files = glob('*.pvd')
 changepvdtimes(pvd_files,t)
 renamevtkfiles(pvd_files)