diff --git a/StructureCode/StructureModel_impl.h b/StructureCode/StructureModel_impl.h
index c8c791f..1210cee 100644
--- a/StructureCode/StructureModel_impl.h
+++ b/StructureCode/StructureModel_impl.h
@@ -1479,7 +1479,26 @@ class StructureModel<T>::Impl
       tractionZ[n][2] = wgPlusTranspose[2][2]*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*eta1[n]
                         -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvalue[n][0]*eigenvector1[n][2]*eigenvector1[n][2] + eigenvalue[n][1]*eigenvector2[n][2]*eigenvector2[n][2] + eigenvalue[n][2]*eigenvector3[n][2]*eigenvector3[n][2]);*/	  
 
-	  tractionX[n][0] = (wgPlusTranspose[0][0]-(wg[0][0]+wg[1][1]+wg[2][2])*2.0/3.0)*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*(eta1[n]+2.0/3.0*eta[n])
+	  tractionX[n][0] = wgPlusTranspose[0][0]*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*eta1[n]
+	                    -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector1[n][0]-(wg[0][0]+wg[1][1]+wg[2][2])/3.0);
+	  tractionX[n][1] = wgPlusTranspose[0][1]*eta[n]
+	                    -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector1[n][1]);
+	  tractionX[n][2] = wgPlusTranspose[0][2]*eta[n]
+	                    -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector1[n][2]);
+      tractionY[n][0] = wgPlusTranspose[1][0]*eta[n]
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector2[n][0]);
+      tractionY[n][1] = wgPlusTranspose[1][1]*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*eta1[n]
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector2[n][1]-(wg[0][0]+wg[1][1]+wg[2][2])/3.0);
+      tractionY[n][2] = wgPlusTranspose[1][2]*eta[n]
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector2[n][2]);
+      tractionZ[n][0] = wgPlusTranspose[2][0]*eta[n]
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector3[n][0]);
+      tractionZ[n][1] = wgPlusTranspose[2][1]*eta[n]
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector3[n][1]);
+      tractionZ[n][2] = wgPlusTranspose[2][2]*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*eta1[n]
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(eigenvector3[n][2]-(wg[0][0]+wg[1][1]+wg[2][2])/3.0);
+
+	  /*tractionX[n][0] = (wgPlusTranspose[0][0]-(wg[0][0]+wg[1][1]+wg[2][2])*2.0/3.0)*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*(eta1[n]+2.0/3.0*eta[n])
 	                    -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(wgPlusTranspose[0][0]/2.0-(wg[0][0]+wg[1][1]+wg[2][2])/3.0);
 	  tractionX[n][1] = wgPlusTranspose[0][1]*eta[n]
 	                    -2.0*etaold[n]*(1-pfv[n]*pfv[n])*wgPlusTranspose[0][1];
@@ -1496,7 +1515,7 @@ class StructureModel<T>::Impl
       tractionZ[n][1] = wgPlusTranspose[2][1]*eta[n]
                         -2.0*etaold[n]*(1-pfv[n]*pfv[n])*wgPlusTranspose[2][1];
       tractionZ[n][2] = (wgPlusTranspose[2][2]-(wg[0][0]+wg[1][1]+wg[2][2])*2.0/3.0)*eta[n]+(wg[0][0]+wg[1][1]+wg[2][2])*(eta1[n]+2.0/3.0*eta[n])
-                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(wgPlusTranspose[2][2]/2.0-(wg[0][0]+wg[1][1]+wg[2][2])/3.0);	  
+                        -2.0*etaold[n]*(1-pfv[n]*pfv[n])*(wgPlusTranspose[2][2]/2.0-(wg[0][0]+wg[1][1]+wg[2][2])/3.0);	  */
 
       
       if ((wg[0][0]+wg[1][1]+wg[2][2])>0 && pfperfect[n]!=-1){
diff --git a/StructureCode/StructureSourceDiscretization.h b/StructureCode/StructureSourceDiscretization.h
index 579a6c5..66927ff 100644
--- a/StructureCode/StructureSourceDiscretization.h
+++ b/StructureCode/StructureSourceDiscretization.h
@@ -315,7 +315,7 @@ class StructureSourceDiscretization : public Discretization
         const T diffMetric = faceAreaMag[f]*faceAreaMag[f]/dot(faceArea[f],ds);
         const VectorT3 secondaryCoeff = faceMu*(faceArea[f]-ds*diffMetric);
         
-        const T divUc0 = (vGradCell[c0][0][0] +vGradCell[c0][1][1] +vGradCell[c0][2][2]);
+        /*const T divUc0 = (vGradCell[c0][0][0] +vGradCell[c0][1][1] +vGradCell[c0][2][2]);
         const T divUc1 = (vGradCell[c1][0][0] +vGradCell[c1][1][1] +vGradCell[c1][2][2]);
         
         faceEigenvalue11[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*(vGradCell[c0][0][0]-divUc0/3.0)*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*(vGradCell[c1][0][0]-divUc1/3.0)*wt1;
@@ -328,7 +328,22 @@ class StructureSourceDiscretization : public Discretization
         
         faceEigenvalue31[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*vGradCell[c0][2][0]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*vGradCell[c1][2][0]*wt1;
         faceEigenvalue32[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*vGradCell[c0][2][1]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*vGradCell[c1][2][1]*wt1;
-        faceEigenvalue33[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*(vGradCell[c0][2][2]-divUc0/3.0)*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*(vGradCell[c1][2][2]-divUc1/3.0)*wt1;
+        faceEigenvalue33[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*(vGradCell[c0][2][2]-divUc0/3.0)*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*(vGradCell[c1][2][2]-divUc1/3.0)*wt1;*/
+
+        const T divUc0 = (eigenvector1Cell[c0][0] +eigenvector2Cell[c0][1] +eigenvector3Cell[c0][2]);
+        const T divUc1 = (eigenvector1Cell[c1][0] +eigenvector2Cell[c1][1] +eigenvector3Cell[c1][2]);
+        
+        faceEigenvalue11[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*(eigenvector1Cell[c0][0]-divUc0/3.0)*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*(eigenvector1Cell[c1][0]-divUc1/3.0)*wt1;
+        faceEigenvalue12[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*eigenvector1Cell[c0][1]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*eigenvector1Cell[c1][1]*wt1;
+        faceEigenvalue13[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*eigenvector1Cell[c0][2]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*eigenvector1Cell[c1][2]*wt1;
+        
+        faceEigenvalue21[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*eigenvector2Cell[c0][0]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*eigenvector2Cell[c1][0]*wt1;
+        faceEigenvalue22[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*(eigenvector2Cell[c0][1]-divUc0/3.0)*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*(eigenvector2Cell[c1][1]-divUc1/3.0)*wt1;
+        faceEigenvalue23[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*eigenvector2Cell[c0][2]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*eigenvector2Cell[c1][2]*wt1;
+        
+        faceEigenvalue31[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*eigenvector3Cell[c0][0]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*eigenvector3Cell[c1][0]*wt1;
+        faceEigenvalue32[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*eigenvector3Cell[c0][1]*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*eigenvector3Cell[c1][1]*wt1;
+        faceEigenvalue33[0]=2.0*(1.0-pfvCell[c0]*pfvCell[c0])*(eigenvector3Cell[c0][2]-divUc0/3.0)*wt0 + 2.0*(1.0-pfvCell[c1]*pfvCell[c1])*(eigenvector3Cell[c1][2]-divUc1/3.0)*wt1;
 
         // mu*grad U ^ T + lambda * div U I
 	source[0] = faceMu*(gradF[0][0]*Af[0] + gradF[0][1]*Af[1] + gradF[0][2]*Af[2])
@@ -352,6 +367,7 @@ class StructureSourceDiscretization : public Discretization
     source[1] -= faceMuOld*(faceEigenvalue32[0]*Af[2] +faceEigenvalue32[1]*Af[2] + faceEigenvalue32[2]*Af[2]);
     source[2] -= faceMuOld*(faceEigenvalue33[0]*Af[2] +faceEigenvalue33[1]*Af[2] + faceEigenvalue33[2]*Af[2]);*/
     
+    
     source[0] -= faceMuOld*(faceEigenvalue11[0]*Af[0] +faceEigenvalue21[0]*Af[1] + faceEigenvalue31[0]*Af[2]);
     source[1] -= faceMuOld*(faceEigenvalue12[0]*Af[0] +faceEigenvalue22[0]*Af[1] + faceEigenvalue32[0]*Af[2]);
     source[2] -= faceMuOld*(faceEigenvalue13[0]*Af[0] +faceEigenvalue23[0]*Af[1] + faceEigenvalue33[0]*Af[2]);
diff --git a/StructureCode/VTKWriter.h b/StructureCode/VTKWriter.h
deleted file mode 100644
index 05242b4..0000000
--- a/StructureCode/VTKWriter.h
+++ /dev/null
@@ -1,436 +0,0 @@
-// This file os part of FVM
-// Copyright (c) 2012 FVM Authors
-// See LICENSE file for terms.
-
-#ifndef _VTKWRITER_H_
-#define _VTKWRITER_H_
-
-#include "atype.h"
-
-#include "Mesh.h"
-#include "ArrayWriter.h"
-#include "Field.h"
-#include "GeomFields.h"
-
-// taken from VTK source
-#define VTK_EMPTY_CELL        0
-#define VTK_VERTEX            1
-#define VTK_POLY_VERTEX       2
-#define VTK_LINE              3
-#define VTK_POLY_LINE         4
-#define VTK_TRIANGLE          5
-#define VTK_TRIANGLE_STRIP    6
-#define VTK_POLYGON           7
-#define VTK_PIXEL             8
-#define VTK_QUAD              9
-#define VTK_TETRA             10
-#define VTK_VOXEL             11
-#define VTK_HEXAHEDRON        12
-#define VTK_WEDGE             13
-#define VTK_PYRAMID           14
-#define VTK_PENTAGONAL_PRISM  15
-#define VTK_HEXAGONAL_PRISM   16
-#define VTK_CONVEX_POINT_SET  41
-
-template<class T>
-class VTKWriter
-{
-public:
-  typedef Array<T> TArray;
-  typedef Array<int> IntArray;
-  typedef Vector<T,3> VectorT3;
-  typedef Array<Vector<T,3> > VectorT3Array;
-  
-  VTKWriter(const GeomFields& geomFields,
-            const MeshList& meshes,
-            const string fileName,
-            const string& comment,
-            const bool binary,
-            const int atypeComponent,
-            const bool surfaceOnly=false) :
-    _geomFields(geomFields),
-    _meshes(meshes),
-    _fp(fopen(fileName.c_str(),"wb")),
-    _binary(binary),
-    _atypeComponent(atypeComponent),
-    _surfaceOnly(surfaceOnly)
-  {
-    if (!_fp)
-      throw CException("VTKWriter: cannot open file " + fileName +
-                       "for writing");
-
-    Field globalIndexField("gIndex");
-
-    const int numMeshes = _meshes.size();
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-        const StorageSite& nodes = mesh.getNodes();
-        const int numNodes = nodes.getCount();
-        shared_ptr<IntArray> gnPtr(new IntArray(numNodes));
-        globalIndexField.addArray(nodes,gnPtr);
-        *gnPtr = -1;
-    }
-
-    _gNodeCount = 0;
-    _gCellCount = 0;
-    int gCellNodeCount =0;
-    
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-        const StorageSite& nodes = mesh.getNodes();
-        IntArray& gNodeIndex = dynamic_cast<IntArray&>(globalIndexField[nodes]);
-
-        if (_surfaceOnly)
-        {
-            foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
-            {
-                const FaceGroup& fg = *fgPtr;
-                if (fg.groupType!="interior")
-                {
-                    const StorageSite& faces = fg.site;
-                    const int faceCount = faces.getCount();
-                    const CRConnectivity& faceNodes = mesh.getFaceNodes(faces);
-
-                    for(int f=0; f<faceCount; f++)
-                    {
-                        const int nFaceNodes = faceNodes.getCount(f);
-                        for(int nn=0; nn<nFaceNodes; nn++)
-                        {
-                            const int node = faceNodes(f,nn);
-                            if (gNodeIndex[node] == -1)
-                              gNodeIndex[node] = _gNodeCount++;
-                        }
-                        _gCellCount++;
-                        gCellNodeCount += nFaceNodes+1;
-                        
-                    }
-                }
-            }
-        }
-        else
-        {
-            const StorageSite& cells = mesh.getCells();
-            const int numCells = cells.getSelfCount();
-            const CRConnectivity& cellNodes = mesh.getCellNodes();
-            const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
-            
-            for(int c=0; c<numCells; c++)
-              if (ibType[c] == Mesh::IBTYPE_FLUID)
-              {
-                  const int nCellNodes = cellNodes.getCount(c);
-                  for(int nn=0; nn<nCellNodes; nn++)
-                  {
-                      const int node = cellNodes(c,nn);
-                      if (gNodeIndex[node] == -1)
-                        gNodeIndex[node] = _gNodeCount++;
-                  }
-                  _gCellCount++;
-                  gCellNodeCount += nCellNodes+1;
-              }
-        }
-    }
-
-    Array<Vector<double,3> > gNodeCoords(_gNodeCount);
-
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-        const StorageSite& nodes = mesh.getNodes();
-        const int numNodes = nodes.getCount();
-
-        const IntArray& gNodeIndex =
-          dynamic_cast<const IntArray&>(globalIndexField[nodes]);
-
-        const VectorT3Array& coords =
-            dynamic_cast<const VectorT3Array&>(geomFields.coordinate[nodes]);
-
-        for(int node=0; node<numNodes; node++)
-        {
-            const int gn = gNodeIndex[node];
-            if (gn != -1)
-            {
-                for(int k=0; k<3; k++)
-                  gNodeCoords[gn][k] =   coords[node][k];
-            }
-        }
-    }
-    
-    fprintf(_fp,"# vtk DataFile Version 2.0\n");
-    fprintf(_fp,"%s\n",comment.c_str());
-    //    if (_surfaceOnly)
-    //   fprintf(_fp,"ASCII\nDATASET POLYDATA\n");
-    //else
-      fprintf(_fp,"ASCII\nDATASET UNSTRUCTURED_GRID\n");
-    
-    fprintf(_fp,"POINTS %d double\n", _gNodeCount);
-
-    for(int node=0; node<_gNodeCount; node++)
-      fprintf(_fp,"%12.5le %12.5le %12.5le\n", gNodeCoords[node][0],
-              gNodeCoords[node][1],gNodeCoords[node][2]);
-
-    //if (_surfaceOnly)
-    //  fprintf(_fp,"LINES");
-    //else
-      fprintf(_fp,"CELLS");
-    
-    fprintf(_fp," %d %d\n", _gCellCount, gCellNodeCount);
-
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-        const StorageSite& nodes = mesh.getNodes();
-        const IntArray& gNodeIndex =
-          dynamic_cast<const IntArray&>(globalIndexField[nodes]);
-
-        if (_surfaceOnly)
-        {
-            foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
-            {
-                const FaceGroup& fg = *fgPtr;
-                if (fg.groupType!="interior")
-                {
-                    const StorageSite& faces = fg.site;
-                    const int faceCount = faces.getCount();
-                    const CRConnectivity& faceNodes = mesh.getFaceNodes(faces);
-
-                    for(int f=0; f<faceCount; f++)
-                    {
-                        const int nFaceNodes = faceNodes.getCount(f);
-                        fprintf(_fp, "%d ", nFaceNodes);
-                        for(int nn=0; nn<nFaceNodes; nn++)
-                        {
-                            const int node = faceNodes(f,nn);
-                            fprintf(_fp, "%d ", gNodeIndex[node]);
-                        }
-                        fprintf(_fp,"\n");
-                    }
-                }
-            }
-        }
-        else
-        {
-            
-            const StorageSite& cells = mesh.getCells();
-            const int numCells = cells.getSelfCount();
-            const IntArray& ibType =
-              dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
-            
-            const CRConnectivity& cellNodes = mesh.getCellNodes();
-            
-            for(int c=0; c<numCells; c++)
-              if (ibType[c] == Mesh::IBTYPE_FLUID)
-              {
-                  const int nCellNodes = cellNodes.getCount(c);
-                  fprintf(_fp, "%d ", nCellNodes);
-                  for(int nn=0; nn<nCellNodes; nn++)
-                  {
-                      const int node = cellNodes(c,nn);
-                      fprintf(_fp, "%d ", gNodeIndex[node]);
-                  }
-                  fprintf(_fp,"\n");
-              }
-        }
-    }
-
-    fprintf(_fp,"CELL_TYPES %d\n", _gCellCount);
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-        const int dim = mesh.getDimension();
-        if (_surfaceOnly)
-        {
-            foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
-            {
-                const FaceGroup& fg = *fgPtr;
-                if (fg.groupType!="interior")
-                {
-                    const StorageSite& faces = fg.site;
-                    const int faceCount = faces.getCount();
-                    const CRConnectivity& faceNodes = mesh.getFaceNodes(faces);
-                    
-                    for(int f=0; f<faceCount; f++)
-                    {
-                        const int nFaceNodes = faceNodes.getCount(f);
-                        int vtkCellType = VTK_POLYGON;
-                        if (dim == 2)
-                          vtkCellType = VTK_LINE;
-                        else if (nFaceNodes == 3)
-                          vtkCellType = VTK_TRIANGLE;
-                        else if (nFaceNodes == 4)
-                          vtkCellType = VTK_QUAD;
-                        fprintf(_fp,"%d\n",vtkCellType);
-                    }
-                }
-            }
-        }
-        else
-        {
-            const StorageSite& cells = mesh.getCells();
-            const int numCells = cells.getSelfCount();
-            const IntArray& ibType =
-              dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
-            
-            const CRConnectivity& cellNodes = mesh.getCellNodes();
-            
-            for(int c=0; c<numCells; c++)
-              if (ibType[c] == Mesh::IBTYPE_FLUID)
-              {
-                  const int nCellNodes = cellNodes.getCount(c);
-                  int vtkCellType = VTK_CONVEX_POINT_SET;
-                  if (dim == 2)
-                  {
-                      if (nCellNodes == 4)
-                        vtkCellType = VTK_QUAD;
-                      else if (nCellNodes == 3)
-                        vtkCellType = VTK_TRIANGLE;
-                      else
-                        vtkCellType = VTK_POLYGON;
-                  }
-                  else
-                  {
-                      if (nCellNodes == 4)
-                        vtkCellType = VTK_TETRA;
-                      else if (nCellNodes == 8)
-                        vtkCellType = VTK_HEXAHEDRON;
-                      else if (nCellNodes == 5)
-                        vtkCellType = VTK_PYRAMID;
-                      else if (nCellNodes == 6)
-                        vtkCellType = VTK_WEDGE;
-                  }
-                  fprintf(_fp,"%d\n",vtkCellType);
-              }
-        }
-    }
-    
-    fprintf(_fp,"CELL_DATA %d\n", _gCellCount);
-
-  }
-
-  void init()
-  {}
-  
-
-  void writeScalarField(const Field& field,
-                        const string label)
-  {
-    const int numMeshes = _meshes.size();
-    fprintf(_fp,"SCALARS %s double\n", label.c_str());
-    fprintf(_fp,"LOOKUP_TABLE default\n");
-    
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-        if (_surfaceOnly)
-        {
-            const StorageSite& allFaces = mesh.getFaces();
-            const TArray& a =
-              dynamic_cast<const TArray&>(field[allFaces]);
-
-            foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
-            {
-                const FaceGroup& fg = *fgPtr;
-                if (fg.groupType!="interior")
-                {
-                    const StorageSite& faces = fg.site;
-                    const int numFaces = faces.getSelfCount();
-                    const int offset = faces.getOffset();
-                    for(int f=0; f<numFaces; f++)
-                    {
-                        const int gf = f + offset;
-                        fprintf(_fp,"%12.5le\n", (a[gf]));
-                    }
-                }
-            }
-        }
-        else
-        {
-            const StorageSite& cells = mesh.getCells();
-            const int numCells = cells.getSelfCount();
-            const IntArray& ibType =
-              dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
-            
-            const TArray& aCell = dynamic_cast<const TArray&>(field[cells]);
-            for(int c=0; c<numCells; c++)
-              if (ibType[c] == Mesh::IBTYPE_FLUID)
-              {
-                  fprintf(_fp,"%12.5le\n", (aCell[c]));
-              }
-        }
-    }
-  }
-  
-  void writeVectorField(const Field& field,const string label)
-  {
-
-
-    fprintf(_fp,"VECTORS %s double\n", label.c_str());
-    const int numMeshes = _meshes.size();
-    
-    for (int n=0; n<numMeshes; n++)
-    {
-        const Mesh& mesh = *_meshes[n];
-
-        if (_surfaceOnly)
-        {
-            const StorageSite& allFaces = mesh.getFaces();
-            const VectorT3Array& v =
-              dynamic_cast<const VectorT3Array&>(field[allFaces]);
-
-            foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
-            {
-                const FaceGroup& fg = *fgPtr;
-                if (fg.groupType!="interior")
-                {
-                    const StorageSite& faces = fg.site;
-                    const int numFaces = faces.getSelfCount();
-                    const int offset = faces.getOffset();
-                    for(int f=0; f<numFaces; f++)
-                    {
-                        const int gf = f + offset;
-                        fprintf(_fp,"%12.5le %12.5le %12.5le\n", (v[gf][0]),
-                                (v[gf][1]),
-                                (v[gf][2])
-                                );
-                    }
-                }
-            }
-        }
-        else
-        {
-            const StorageSite& cells = mesh.getCells();
-            const int numCells = cells.getSelfCount();
-            const IntArray& ibType =
-              dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
-            
-            const VectorT3Array& aCell =
-              dynamic_cast<const VectorT3Array&>(field[cells]);
-            for(int c=0; c<numCells; c++)
-              if (ibType[c] == Mesh::IBTYPE_FLUID)
-              {
-                  fprintf(_fp,"%12.5le %12.5le %12.5le\n", (aCell[c][0]),
-                          (aCell[c][1]),
-                          (aCell[c][2])
-                          );
-              }
-        }
-    }
-  }
-
-  void finish()
-  {
-    fclose(_fp);
-  }
-
-private:
-  const GeomFields& _geomFields;
-  const MeshList _meshes;
-  FILE *_fp;
-  const bool _binary;
-  const int _atypeComponent;
-  int _gNodeCount;
-  int _gCellCount;
-  const bool _surfaceOnly;
-};
-#endif
diff --git a/StructureCode/VTKWriter.i b/StructureCode/VTKWriter.i
deleted file mode 100644
index ec1143f..0000000
--- a/StructureCode/VTKWriter.i
+++ /dev/null
@@ -1,35 +0,0 @@
-%{
-#include "VTKWriter.h"
-#include "atype.h"
-%}
-
-%include "std_string.i"
-
-%include "atype.i"
-%import "Mesh.i"
-%import "GeomFields.h"
-
-using namespace std;
-
-template<class T>
-class VTKWriter
-{
-public:
-  VTKWriter(const GeomFields& geomFields,
-            const MeshList& meshes,
-            const string fileName,
-            const string comment,
-            const bool binary,
-            const int atypeComponent,
-            const bool surfaceOnly=false);
-  void init();
-  void finish();
-  void writeScalarField(const Field& field,
-                        const string label);
-  void writeVectorField(const Field& field,
-                        const string label);
-};
-
-%template(VTKWriterA) VTKWriter< ATYPE_STR >;
-
-
diff --git a/fixedFiber/fixedFiber.py b/fixedFiber/fixedFiber.py
index 3f7b394..91d55c8 100644
--- a/fixedFiber/fixedFiber.py
+++ b/fixedFiber/fixedFiber.py
@@ -753,7 +753,11 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
                    
                    #decomposeStrainTensor(strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i],eigenvalueFieldsA[i],eigenvector1FieldsA[i],eigenvector2FieldsA[i],eigenvector3FieldsA[i],\
                    #i,pfperfectFieldsA[i],rank_id)
-           
+                   for j in range (0,3):
+                       eigenvector1FieldsA[i][j]=strainXFieldsA[i][j]
+                       eigenvector2FieldsA[i][j]=strainYFieldsA[i][j]
+                       eigenvector3FieldsA[i][j]=strainZFieldsA[i][j]
+                   
                    if strain_trace[i] >= 0 and SymFlag==2:
                        if V_flag[i] == 1 :
                            #if strain_trace[i] < 1e-1 :
@@ -878,6 +882,9 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            strain_dev2_trace=(strainXFieldsA[i][0]-strain_trace_mean)**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
            strainYFieldsA[i][0]**2+(strainYFieldsA[i][1]-strain_trace_mean)**2+strainYFieldsA[i][2]**2+\
            strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2]-strain_trace_mean)**2
+           strain_2_trace=(strainXFieldsA[i][0])**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
+           strainYFieldsA[i][0]**2+(strainYFieldsA[i][1])**2+strainYFieldsA[i][2]**2+\
+           strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2])**2
                #if (K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace) > EnergyHistoryField[i]:
                #    ElasticEnergyField[i]=K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
                #    #EnergyHistoryField[i]=ElasticEnergyField[i]
@@ -906,8 +913,12 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            else:
                if strain_trace[i] >0:
                    ElasticEnergyField[i] = K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_2_trace
                else: 
                    ElasticEnergyField[i] = G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = G_local[i]*strain_2_trace
                #if i==100:
                #    print i,ElasticEnergyField[i],eigenvalueFieldsA[i]
                #    print strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i]
diff --git a/homoCase/homoCase.py b/homoCase/homoCase.py
index d8d5c5a..f6283b4 100644
--- a/homoCase/homoCase.py
+++ b/homoCase/homoCase.py
@@ -226,7 +226,7 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
 DeformUnit = (cFED*BoundaryPositionTop/Lamda)**0.5  #Normalized Displacement Unit
 #DispStep = 0.02*DeformUnit	   # Displacement Step
 DispStep = 1e-6
-StressStep = 1e6
+StressStep = 2e6
 LoadCoef = -0.5
 
 OInterval_s = 1                 #Output interval for equilibrium status
@@ -239,8 +239,8 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
 MidIterUpLimit = 200
 
 StiffnessResidual = 1e-6       #Used to have a lower bound of the material constant for damaged cell
-StructTolerance = 1e-3         #Tolerance for structure model inner iteration
-StructOuterTolerance = 1e-3
+StructTolerance = 1e-4         #Tolerance for structure model inner iteration
+StructOuterTolerance = 1e-4
 StructIterFlag = 0             #1--Do structure model iteration; 0--No structure model iteration
 StructIterUpLimit = 80
 
@@ -765,7 +765,11 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
                    
                    #decomposeStrainTensor(strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i],eigenvalueFieldsA[i],eigenvector1FieldsA[i],eigenvector2FieldsA[i],eigenvector3FieldsA[i],\
                    #i,pfperfectFieldsA[i],rank_id)
-           
+                   for j in range (0,3):
+                       eigenvector1FieldsA[i][j]=strainXFieldsA[i][j]
+                       eigenvector2FieldsA[i][j]=strainYFieldsA[i][j]
+                       eigenvector3FieldsA[i][j]=strainZFieldsA[i][j]
+                   
                    if strain_trace[i] >= 0 and SymFlag==2:
                        if V_flag[i] == 1 :
                            #if strain_trace[i] < 1e-1 :
@@ -887,6 +891,9 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            strain_dev2_trace=(strainXFieldsA[i][0]-strain_trace_mean)**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
            strainYFieldsA[i][0]**2+(strainYFieldsA[i][1]-strain_trace_mean)**2+strainYFieldsA[i][2]**2+\
            strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2]-strain_trace_mean)**2
+           strain_2_trace=(strainXFieldsA[i][0])**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
+           strainYFieldsA[i][0]**2+(strainYFieldsA[i][1])**2+strainYFieldsA[i][2]**2+\
+           strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2])**2
                #if (K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace) > EnergyHistoryField[i]:
                #    ElasticEnergyField[i]=K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
                #    #EnergyHistoryField[i]=ElasticEnergyField[i]
@@ -915,8 +922,12 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            else:
                if strain_trace[i] >0:
                    ElasticEnergyField[i] = K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_2_trace
                else: 
                    ElasticEnergyField[i] = G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = G_local[i]*strain_2_trace
                if i==100:
                    print i,ElasticEnergyField[i],eigenvalueFieldsA[i]
                    print strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i]
diff --git a/multiFiber-3D/multiFiber.py b/multiFiber-3D/multiFiber.py
index 4dd1931..8640195 100644
--- a/multiFiber-3D/multiFiber.py
+++ b/multiFiber-3D/multiFiber.py
@@ -229,24 +229,24 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
 DispStep = 0.03*DeformUnit	   # Displacement Step
 StressStep = 1e6
 
-OInterval_s = 1                 #Output interval for equilibrium status
+OInterval_s = 4                 #Output interval for equilibrium status
 OInterval_l = 20
 MidOInterval_s = 50              #Output interval for intermediate status
 MidOInterval_l = 50
-OPFLimit = 0.02
+OPFLimit = 0.04
 OUpLimit=40                   #Upper Limit for large displacement step
 DispReFactor=1.0               #Smaller displacement step is: 1/DispReFactor of larger displacement step
 MidIterUpLimit = 80
 
 StiffnessResidual = 1e-5       #Used to have a lower bound of the material constant for damaged cell
-StructTolerance = 5e-3         #Tolerance for structure model inner iteration
-StructOuterTolerance = 5e-3
+StructTolerance = 1e-3         #Tolerance for structure model inner iteration
+StructOuterTolerance = 1e-3
 StructIterFlag = 1             #1--Do structure model iteration; 0--No structure model iteration
 StructIterUpLimit = 50
 
-PFTolerance = 5e-3             #Tolerance for fracture model iteration
-PFOuterTolerance = 5e-3
-PFIterFlag = 0                 #1--Do convergence test iteration; 0--No convergence test iteration
+PFTolerance = 1e-3             #Tolerance for fracture model iteration
+PFOuterTolerance = 1e-3
+PFIterFlag = 1                 #1--Do convergence test iteration; 0--No convergence test iteration
 
 PerfectRad = 0e-3
 SymFlag = 0 # 1--Symmetric 0--Asymmetric
@@ -823,7 +823,11 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
                    
                    #decomposeStrainTensor(strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i],eigenvalueFieldsA[i],eigenvector1FieldsA[i],eigenvector2FieldsA[i],eigenvector3FieldsA[i],\
                    #i,pfperfectFieldsA[i],rank_id)
-           
+                   for j in range (0,3):
+                       eigenvector1FieldsA[i][j]=strainXFieldsA[i][j]
+                       eigenvector2FieldsA[i][j]=strainYFieldsA[i][j]
+                       eigenvector3FieldsA[i][j]=strainZFieldsA[i][j]
+                   
                    if strain_trace[i] >= 0 and SymFlag==2:
                        if V_flag[i] == 1 :
                            #if strain_trace[i] < 1e-1 :
@@ -948,6 +952,9 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            strain_dev2_trace=(strainXFieldsA[i][0]-strain_trace_mean)**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
            strainYFieldsA[i][0]**2+(strainYFieldsA[i][1]-strain_trace_mean)**2+strainYFieldsA[i][2]**2+\
            strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2]-strain_trace_mean)**2
+           strain_2_trace=(strainXFieldsA[i][0])**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
+           strainYFieldsA[i][0]**2+(strainYFieldsA[i][1])**2+strainYFieldsA[i][2]**2+\
+           strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2])**2
                #if (K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace) > EnergyHistoryField[i]:
                #    ElasticEnergyField[i]=K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
                #    #EnergyHistoryField[i]=ElasticEnergyField[i]
@@ -976,8 +983,12 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            else:
                if strain_trace[i] >0:
                    ElasticEnergyField[i] = K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_2_trace
                else: 
                    ElasticEnergyField[i] = G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = G_local[i]*strain_2_trace
                #if i==100:
                #    print i,ElasticEnergyField[i],eigenvalueFieldsA[i]
                #    print strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i]
diff --git a/multiFiber-3D/multifiber-3D-coarse.cas b/multiFiber-3D/multifiber-3D-coarse.cas
old mode 100755
new mode 100644
diff --git a/multiFiber/fixedFiber.py b/multiFiber/fixedFiber.py
index e11d3ba..9661db1 100644
--- a/multiFiber/fixedFiber.py
+++ b/multiFiber/fixedFiber.py
@@ -20,8 +20,8 @@
 from optparse import OptionParser
 from FluentCase import FluentCase
 
-import tecplotEntireStructureDomainPara
-import tecplotEntireFractureDomainPara
+import vtkEntireStructureDomainPara
+import vtkEntireFractureDomainPara
 
 def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pfp_flag,rank):
     zeroThreshold1=1e-12
@@ -741,7 +741,11 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
                    
                    #decomposeStrainTensor(strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i],eigenvalueFieldsA[i],eigenvector1FieldsA[i],eigenvector2FieldsA[i],eigenvector3FieldsA[i],\
                    #i,pfperfectFieldsA[i],rank_id)
-           
+                   for j in range (0,3):
+                       eigenvector1FieldsA[i][j]=strainXFieldsA[i][j]
+                       eigenvector2FieldsA[i][j]=strainYFieldsA[i][j]
+                       eigenvector3FieldsA[i][j]=strainZFieldsA[i][j]
+                   
                    if strain_trace[i] >= 0 and SymFlag==2:
                        if V_flag[i] == 1 :
                            #if strain_trace[i] < 1e-1 :
@@ -863,6 +867,9 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            strain_dev2_trace=(strainXFieldsA[i][0]-strain_trace_mean)**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
            strainYFieldsA[i][0]**2+(strainYFieldsA[i][1]-strain_trace_mean)**2+strainYFieldsA[i][2]**2+\
            strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2]-strain_trace_mean)**2
+           strain_2_trace=(strainXFieldsA[i][0])**2+strainXFieldsA[i][1]**2+strainXFieldsA[i][2]**2+\
+           strainYFieldsA[i][0]**2+(strainYFieldsA[i][1])**2+strainYFieldsA[i][2]**2+\
+           strainZFieldsA[i][0]**2+strainZFieldsA[i][1]**2+(strainZFieldsA[i][2])**2
                #if (K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace) > EnergyHistoryField[i]:
                #    ElasticEnergyField[i]=K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
                #    #EnergyHistoryField[i]=ElasticEnergyField[i]
@@ -891,13 +898,17 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
            else:
                if strain_trace[i] >0:
                    ElasticEnergyField[i] = K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = Lamda_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_2_trace
                else: 
                    ElasticEnergyField[i] = G_local[i]*strain_dev2_trace
-               if i==100:
-                   print i,ElasticEnergyField[i],eigenvalueFieldsA[i]
-                   print strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i]
-                   print tractZFieldsA[i][2]
-                   print eigenvector1FieldsA[i],eigenvector2FieldsA[i],eigenvector3FieldsA[i]
+                   #ElasticEnergyField[i] = G_local[i]*(eigenvalueFieldsA[i][0]**2.0+eigenvalueFieldsA[i][1]**2.0+eigenvalueFieldsA[i][2]**2.0)
+                   #ElasticEnergyField[i] = G_local[i]*strain_2_trace
+               #if i==100:
+               #    print i,ElasticEnergyField[i],eigenvalueFieldsA[i]
+               #    print strainXFieldsA[i],strainYFieldsA[i],strainZFieldsA[i]
+               #    print tractZFieldsA[i][2]
+               #    print eigenvector1FieldsA[i],eigenvector2FieldsA[i],eigenvector3FieldsA[i]
                #    #print tractXFieldsA[i][0],tractYFieldsA[i][1],tractZFieldsA[i][2],pfvFieldsA[i],V_flag[i]
                Total_Elastic_Energy[0] = Total_Elastic_Energy[0] + ((PhaseFieldA[i]**2.0+StiffnessResidual)*(K_local[i]/2.0*strain_trace_positive**2+G_local[i]*strain_dev2_trace)+K_local[i]/2.0*strain_trace_negative**2)*volumeA[i]
            
@@ -1111,28 +1122,11 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
        else :
            MidOInterval=MidOInterval_l
        #Output Structure Module 
-       if mid_iter % MidOInterval ==0:
-           tecplotEntireStructureDomainPara.dumpTecplotEntireStructureDomain(nmesh,  meshes, fluent_meshes, options.type, structureFields,structure_file_name,title_name,Total_count) 
+       #if mid_iter % MidOInterval ==0:
        #Output Fracture Module 
-       if mid_iter % MidOInterval ==0:
-           tecplotEntireFractureDomainPara.dumpTecplotEntireFractureDomain(nmesh,  meshes, fluent_meshes, options.type, fractureFields,fracture_file_name,title_name,Total_count)
+       #if mid_iter % MidOInterval ==0:
        #Output equilibrium status VTK files
        if mid_iter % MidOInterval ==0:   
-           writer_fracture_iter = exporters.VTKWriterA(geomFields,meshes,"fracture-inter-"+str(nstep)+"-"+str(mid_iter)+"-"+".vtk","fracture output",False,0)
-           writer_fracture_iter.init()
-           writer_fracture_iter.writeScalarField(fractureFields.phasefieldvalue,"PhaseField")
-           writer_fracture_iter.writeVectorField(structureFields.deformation,"Deformation")
-           writer_fracture_iter.finish()
-           writer_structure_iter = exporters.VTKWriterA(geomFields,meshes,"structure-inter-"+str(nstep)+"-"+str(mid_iter)+"-"+".vtk","structure output",False,0)
-           writer_structure_iter.init()
-           writer_structure_iter.writeVectorField(structureFields.deformation,"Deformation")
-           writer_structure_iter.writeVectorField(structureFields.strainX,"StrainX")
-           writer_structure_iter.writeVectorField(structureFields.strainY,"StrainY")
-           writer_structure_iter.writeVectorField(structureFields.strainZ,"StrainZ")
-           writer_structure_iter.writeVectorField(structureFields.tractionX,"TractionX")
-           writer_structure_iter.writeVectorField(structureFields.tractionY,"TractionY")
-           writer_structure_iter.writeVectorField(structureFields.tractionZ,"TractionZ")       
-           writer_structure_iter.finish() 
            Total_count = Total_count +1
        if rank_id == 0 :
            t1 = time.time()
@@ -1186,28 +1180,12 @@ def decomposeStrainTensor (strX,strY,strZ,evalue,evector1,evector2,evector3,i,pf
        "\n")
        sp_structure.flush()         
    title_name="Equil "+str(nstep)
+   #if nstep % OInterval ==0:   
    #Output Structure Module 
    if nstep % OInterval ==0:
-       tecplotEntireStructureDomainPara.dumpTecplotEntireStructureDomain(nmesh,  meshes, fluent_meshes, options.type, structureFields,structure_file_name,title_name,Total_count) 
+       vtkEntireStructureDomainPara.dumpvtkEntireStructureDomain(geomFields, nmesh,  meshes, fluent_meshes, options.type, structureFields,structure_file_name,title_name,nstep) 
    #Output Fracture Module 
    if nstep % OInterval ==0:
-       tecplotEntireFractureDomainPara.dumpTecplotEntireFractureDomain(nmesh,  meshes, fluent_meshes, options.type, fractureFields,fracture_file_name,title_name,Total_count)
-   #Output equilibrium status VTK files
-   if nstep % OInterval ==0:   
-       writer_fracture = exporters.VTKWriterA(geomFields,meshes,"fracture-equil-"+str(nstep)+".vtk","fracture output",False,0)
-       writer_fracture.init()
-       writer_fracture.writeScalarField(fractureFields.phasefieldvalue,"PhaseField")
-       writer_fracture.writeVectorField(structureFields.deformation,"Deformation")
-       writer_fracture.finish()
-       writer_structure = exporters.VTKWriterA(geomFields,meshes,"structure-equil-"+str(nstep)+".vtk","structure output",False,0)
-       writer_structure.init()
-       writer_structure.writeVectorField(structureFields.deformation,"Deformation")
-       writer_structure.writeVectorField(structureFields.strainX,"StrainX")
-       writer_structure.writeVectorField(structureFields.strainY,"StrainY")
-       writer_structure.writeVectorField(structureFields.strainZ,"StrainZ")
-       writer_structure.writeVectorField(structureFields.tractionX,"TractionX")
-       writer_structure.writeVectorField(structureFields.tractionY,"TractionY")
-       writer_structure.writeVectorField(structureFields.tractionZ,"TractionZ")      
-       writer_structure.finish()  
-       Total_count = Total_count +1
+       vtkEntireFractureDomainPara.dumpvtkEntireFractureDomain(geomFields, nmesh,  meshes, fluent_meshes, options.type, fractureFields, structureFields,fracture_file_name,title_name,nstep)
+
    #End of Output Equilibrium Status