Merge pull request #243 from joakim-hove/init-check

Init check

.gitignore

@@ -29,3 +29,6 @@ test_vec
 
 # Mac OS X debug info
 *.dSYM
+
+# emacs directory setting:
+.dir-locals.el

CMakeLists_files.cmake

@@ -77,6 +77,7 @@ list (APPEND EXAMPLE_SOURCE_FILES
 	examples/sim_2p_comp_ad.cpp
 	examples/sim_2p_incomp_ad.cpp
 	examples/sim_simple.cpp
+        examples/opm_init_check.cpp     
 	)
 
 # programs listed here will not only be compiled, but also marked for
@@ -84,6 +85,7 @@ list (APPEND EXAMPLE_SOURCE_FILES
 list (APPEND PROGRAM_SOURCE_FILES
 	examples/sim_2p_incomp_ad.cpp
 	examples/sim_fibo_ad.cpp
+        examples/opm_init_check.cpp
 	)
 
 # originally generated with the command:

examples/opm_init_check.cpp

@@ -0,0 +1,362 @@
+/*
+  Copyright 2014 Statoil ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <iostream>
+#include <fstream>
+#include <memory>
+#include <vector>
+#include <string>
+
+#include <ert/ecl/ecl_file.h>
+#include <ert/ecl/ecl_kw.h>
+#include <ert/ecl/ecl_grid.h>
+#include <ert/ecl/ecl_nnc_export.h>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+
+
+#include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+
+
+using namespace Opm;
+
+
+
+
+class CellTrans {
+public: 
+    CellTrans(const std::tuple<int,int,int>& ijk) :
+        m_ijk(ijk) 
+    {
+        
+    }
+    
+    void update(const std::tuple<int,int,int>& ijk , double trans) {
+        auto iter = m_trans.find( ijk );
+        if (iter == m_trans.end())
+            m_trans.insert( std::pair<std::tuple<int,int,int> , double>(ijk , trans));
+        else
+            iter->second *= trans;
+    }
+    
+    int numConnections() const {
+        return m_trans.size();
+    }
+    
+
+    static void jsonDumpIJK(std::ostream& os , const std::tuple<int,int,int>& ijk ) {
+        os << "[" << std::get<0>(ijk) << "," << std::get<1>(ijk) << "," << std::get<2>(ijk) << "]";
+    }
+
+    
+    void jsonDump(std::ostream& os , bool first) {
+        if (!first)
+            os <<",";
+
+        
+        os << "[";
+        jsonDumpIJK(os , m_ijk );
+        os << ",[";
+        {
+            size_t count = 0;
+            for (auto transIter = m_trans.begin(); transIter != m_trans.end(); ++transIter) {
+                std::tuple<int,int,int> ijk = transIter->first;
+                double t = transIter->second;
+            
+                os << "[";
+                jsonDumpIJK( os ,  ijk );
+                os << "," << t << "]";
+                count++;
+                
+                if (count < m_trans.size())
+                    os << ",";
+                else
+                    os << "]";
+            }
+        }
+        os << "]" << std::endl;
+    }
+
+
+    std::tuple<int,int,int> m_ijk;
+    std::map<std::tuple<int,int,int> , double> m_trans;
+};
+    
+
+
+
+class TransGraph {
+public:
+    TransGraph(int nx , int ny , int nz) :
+        m_nx(nx),
+        m_ny(ny),
+        m_nz(nz) 
+    {
+        m_transVector.resize( nx*ny*nz );
+    }
+
+
+    void update(int global_index1 , int global_index2 , double trans) {
+        int size = m_nx * m_ny * m_nz;
+        if ((global_index1 >= 0)   && 
+            (global_index2 >= 0)   &&
+            (global_index1 < size) && 
+            (global_index2 < size)) {
+            
+            size_t g1 = std::min( global_index1 , global_index2 );
+            size_t g2 = std::max( global_index1 , global_index2 );
+
+            std::shared_ptr<CellTrans> cellTrans = m_transVector[g1];
+            if (!cellTrans) {
+                cellTrans = std::make_shared<CellTrans>( getIJK(g1) );
+                m_transVector[g1] = cellTrans;
+            }
+            cellTrans->update( getIJK(g2) , trans );
+            
+        }
+    }
+
+
+    int activeCells() const {
+        int count = 0;
+        for (size_t g= 0; g < m_transVector.size(); g++) {
+            std::shared_ptr<CellTrans> cellTrans = m_transVector[g];
+            if (cellTrans)
+                count++;
+        }
+        return count;
+    }
+
+
+    int activeConnections() const {
+        int count = 0;
+        for (size_t g= 0; g < m_transVector.size(); g++) {
+            std::shared_ptr<CellTrans> cellTrans = m_transVector[g];
+            if (cellTrans)
+                count += cellTrans->numConnections();
+        }
+        return count;
+    }
+
+    
+
+    std::tuple<int,int,int> getIJK(int g) const {
+        int k = g / (m_nx * m_ny);
+        int j = (g - k*m_nx*m_ny) / m_nx;
+        int i = g - k*m_nx*m_ny - j*m_nx;
+
+        return std::tuple<int,int,int>(i,j,k);
+    }
+
+
+    void jsonDump(const std::string& outputFile) {
+        std::ofstream os;
+        bool first = true;
+        os.open(outputFile.c_str());
+        
+        os << "{ \"dims\" : [" << m_nx << "," << m_ny << "," << m_nz << "]" << " , \"graph\":[";
+        for (size_t g= 0; g < m_transVector.size(); g++) {
+            std::shared_ptr<CellTrans> cellTrans = m_transVector[g];
+            if (cellTrans) {
+                cellTrans->jsonDump(os , first);
+                first = false;
+            }
+        }
+        os << "]}" << std::endl;
+        os.close();
+    }
+
+
+
+    int m_nx;
+    int m_ny;
+    int m_nz;
+    std::vector<std::shared_ptr<CellTrans> > m_transVector;
+};
+
+
+/*****************************************************************/
+
+void initOPMTrans(TransGraph& opmTrans , DeckConstPtr deck , std::shared_ptr<const EclipseState> eclipseState) {
+    std::shared_ptr<GridManager> grid = std::make_shared<GridManager>( eclipseState->getEclipseGrid(), eclipseState->getDoubleGridProperty("PORV")->getData());
+    const struct UnstructuredGrid * cGrid = grid->c_grid();
+    std::shared_ptr<BlackoilPropsAdInterface> props;
+
+    props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
+    DerivedGeology geology(*grid->c_grid() , *props, eclipseState);
+    const double * opm_trans_data = geology.transmissibility().data();
+    double SIconversion = Opm::unit::cubic(Opm::unit::meter) * Opm::unit::day * Opm::unit::barsa / (Opm::prefix::centi * Opm::unit::Poise);
+    
+    {
+        for (int face_index = 0; face_index < cGrid->number_of_faces; face_index++ ) {
+            int global_index1 = cGrid->global_cell[ cGrid->face_cells[2*face_index] ];
+            int global_index2 = cGrid->global_cell[ cGrid->face_cells[2*face_index + 1] ];
+
+            opmTrans.update( global_index1 , global_index2 , opm_trans_data[ face_index ] * SIconversion );
+        }
+    }
+}
+
+
+void initEclipseTrans(TransGraph& eclipseTrans , const ecl_grid_type * ecl_grid , const ecl_file_type * ecl_init) {
+    int nx = ecl_grid_get_nx( ecl_grid );
+    int ny = ecl_grid_get_ny( ecl_grid );
+    int nz = ecl_grid_get_nz( ecl_grid );
+
+    if (ecl_file_has_kw( ecl_init , "TRANX")) {
+        ecl_kw_type * tranx_kw = ecl_file_iget_named_kw( ecl_init , "TRANX" , 0 );
+        ecl_kw_type * trany_kw = ecl_file_iget_named_kw( ecl_init , "TRANY" , 0 );
+        ecl_kw_type * tranz_kw = ecl_file_iget_named_kw( ecl_init , "TRANZ" , 0 );
+        for (int k=0; k < nz; k++) {
+            for (int j= 0; j < ny; j++) {
+                for (int i=0; i < nx; i++) {
+                    if (ecl_grid_cell_active3( ecl_grid , i , j , k )) {
+                        size_t g1 = ecl_grid_get_global_index3( ecl_grid , i , j , k );
+                        int a = ecl_grid_get_active_index1( ecl_grid , g1 );
+                        if (a >= 0) {
+                            if (i < (nx - 1) && ecl_grid_cell_active3( ecl_grid , i + 1 , j , k)) {
+                                size_t g2 = ecl_grid_get_global_index3( ecl_grid , i + 1, j , k );
+                                eclipseTrans.update( g1 , g2     , ecl_kw_iget_float( tranx_kw , a ));                            
+                            }
+                            
+                            
+                            if (j < (ny - 1) && ecl_grid_cell_active3( ecl_grid , i , j + 1, k)) {
+                                size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j + 1, k );
+                                eclipseTrans.update( g1 , g2     , ecl_kw_iget_float( trany_kw , a ));                            
+                            }
+                            
+                            
+                            if (k < (nz - 1) && ecl_grid_cell_active3( ecl_grid , i , j , k + 1)) {
+                                size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j , k + 1 );
+                                eclipseTrans.update( g1 , g2     , ecl_kw_iget_float( tranz_kw , a ));                            
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    } else
+        std::cerr << "Init file does not have TRAN[XYZ] keywords" << std::endl;
+    
+    if (ecl_file_has_kw( ecl_init , "TRANX-")) {
+        ecl_kw_type * tranxm_kw = ecl_file_iget_named_kw( ecl_init , "TRANX-" , 0 );
+        ecl_kw_type * tranym_kw = ecl_file_iget_named_kw( ecl_init , "TRANY-" , 0 );
+        ecl_kw_type * tranzm_kw = ecl_file_iget_named_kw( ecl_init , "TRANZ-" , 0 );
+        for (int k=0; k < nz; k++) {
+            for (int j= 0; j < ny; j++) {
+                for (int i=0; i < nx; i++) {
+                    if (ecl_grid_cell_active3( ecl_grid , i , j , k )) {
+                        size_t g1 = ecl_grid_get_global_index3( ecl_grid , i , j , k );
+                        int a = ecl_grid_get_active_index1( ecl_grid , g1 );
+
+                        if (a >= 0) {
+                            if (i > 0 && ecl_grid_cell_active3( ecl_grid , i - 1 , j , k)) {
+                                size_t g2 = ecl_grid_get_global_index3( ecl_grid , i - 1, j , k );
+                                eclipseTrans.update( g1 , g2     , ecl_kw_iget_float( tranxm_kw , a ));                            
+                            }
+                            
+                            
+                            if (j > 0 && ecl_grid_cell_active3( ecl_grid , i , j - 1, k)) {
+                                size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j - 1, k );
+                                eclipseTrans.update( g1 , g2     , ecl_kw_iget_float( tranym_kw , a ));                            
+                            }
+                            
+                            
+                            if (k > 0 && ecl_grid_cell_active3( ecl_grid , i , j , k - 1)) {
+                                size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j , k - 1 );
+                                eclipseTrans.update( g1 , g2     , ecl_kw_iget_float( tranzm_kw , a ));                            
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // NNC
+    {
+        size_t num_nnc = static_cast<size_t>( ecl_nnc_export_get_size( ecl_grid ));
+        std::vector<ecl_nnc_type> nnc(num_nnc);
+        
+        ecl_nnc_export( ecl_grid , ecl_init , nnc.data());
+        for (auto nnc_iter = nnc.begin(); nnc_iter != nnc.end(); ++nnc_iter) 
+            eclipseTrans.update( nnc_iter->global_index1 , nnc_iter->global_index2 , nnc_iter->trans );
+    }
+}
+
+
+
+void dump_transGraph( DeckConstPtr deck , std::shared_ptr<const EclipseState> eclipseState , const ecl_grid_type * ecl_grid , const ecl_file_type * ecl_init , size_t verbosity) {
+    int nx = ecl_grid_get_nx( ecl_grid );
+    int ny = ecl_grid_get_ny( ecl_grid );
+    int nz = ecl_grid_get_nz( ecl_grid );
+    TransGraph opmTrans(nx , ny , nz );
+    TransGraph eclipseTrans( nx , ny , nz);
+    
+    initOPMTrans( opmTrans , deck , eclipseState );
+    initEclipseTrans( eclipseTrans , ecl_grid , ecl_init );
+    opmTrans.jsonDump("opm_trans.json");
+    eclipseTrans.jsonDump("eclipse_trans.json");
+}
+
+
+
+int main(int argc, char** argv) {
+    if (argc < 4) {
+        std::cerr << "The opm_init_check program needs three arguments:" << std::endl <<  std::endl;;
+        std::cerr << "    ECLIPSE.DATA  ECLIPSE.INIT   ECLIPSE.EGRID" << std::endl << std::endl;
+        std::cerr << "Where the ECLIPSE.INIT and ECLIPSE.EGRID are existing binary files";
+        exit(1);
+    }
+
+    std::string input_file = argv[1];
+    std::string init_file = argv[2];
+    std::string grid_file = argv[3];
+    
+    ParserPtr parser(new Parser());
+
+    std::cout << "Parsing input file ............: " << input_file << std::endl;
+    DeckConstPtr deck = parser->parseFile(input_file, false);
+    std::shared_ptr<EclipseState> state = std::make_shared<EclipseState>( deck );
+    
+    std::cout << "Loading eclipse INIT file .....: " << init_file << std::endl;
+    ecl_file_type * ecl_init = ecl_file_open( init_file.c_str() , 0 );
+
+    std::cout << "Loading eclipse EGRID file ....: " << grid_file << std::endl;
+    ecl_grid_type * ecl_grid = ecl_grid_alloc( grid_file.c_str() );
+
+    dump_transGraph( deck , state , ecl_grid , ecl_init , 3);
+    
+    ecl_file_close( ecl_init );
+    ecl_grid_free( ecl_grid );
+    return 0;
+}
+

examples/python/example.py

@@ -0,0 +1,136 @@
+#!/usr/bin/env python
+import sys
+from trans_graph import TransGraph
+
+# This file just contains some example functions for how one can poke
+# around in the TransGraph datastructure.
+
+
+def direction_count(tg):
+    dir_count = {"X" : 0 , 
+                 "Y" : 0 ,
+                 "Z" : 0 ,
+                 "NNC" : 0 }
+    
+    for cell in tg:
+        if cell:
+            for conn in cell:
+                dir_count[ conn.dir ] += 1
+    
+    dir_count["Total"] = dir_count["X"] + dir_count["Y"] + dir_count["Z"] + dir_count["NNC"]
+    return dir_count
+                
+
+
+
+def print_cell( prefix , cell ):
+    print "%s: Cell: (%d,%d,%d) " % (prefix , cell.i , cell.j , cell.k)
+    for conn in cell:
+        print "    Connection => (%3d,%3d,%3d)   Transmissibility: %g     Direction: %s" % (conn.i , conn.j , conn.k , conn.T , conn.dir)
+    
+    
+    print "    cell[\"X\"] => %s" % cell["X"]
+    print "    cell[\"Y\"] => %s" % cell["Y"]
+    print "    cell[\"Z\"] => %s" % cell["Z"]
+    
+
+def connection_to_count(tg , i,j,k):
+    count = 0
+    for cell in tg:
+        if cell.connectsWith(i,j,k):
+            count += 1
+    
+    return count
+        
+
+
+
+
+
+
+
+
+
+
+
+#-----------------------------------------------------------------
+
+def direction_example( opm_tg , ecl_tg ):
+    opm_count = direction_count( opm_tg )
+    ecl_count = direction_count( ecl_tg )
+
+    print "OPM: %s" % opm_count
+    print "ECL: %s" % ecl_count
+    print
+
+
+
+def cell_example(opm_tg , ecl_tg ):
+    opm_cell = opm_tg[21,27,10]
+    ecl_cell = ecl_tg[21,27,10]
+
+    print_cell( "OPM: " , opm_cell )
+    print_cell( "ECL: " , ecl_cell )
+
+        
+
+
+def count_example(opm_tg , ecl_tg ):
+    i = 5
+    j = 10
+    k = 7
+    print "Opm connections to (%d,%d,%d): %d" % (i,j,k , connection_to_count(opm_tg , i,j,k))
+    print "Ecl connections to (%d,%d,%d): %d" % (i,j,k , connection_to_count(ecl_tg , i,j,k))
+    print 
+
+
+def xtrace_example( opm_tg , ecl_tg , j , k):
+    opm_trace = opm_tg.getXTrace( j , k)
+    ecl_trace = ecl_tg.getXTrace( j , k)
+
+    print "OPM: %s" % opm_trace
+    print "ECL: %s" % ecl_trace
+    print
+
+
+def ytrace_example( opm_tg , ecl_tg , i , k):
+    opm_trace = opm_tg.getYTrace(  i , k )
+    ecl_trace = ecl_tg.getYTrace(  i , k )
+
+    print "OPM: %s" % opm_trace
+    print "ECL: %s" % ecl_trace
+    print
+
+
+def ztrace_example( opm_tg , ecl_tg , i , j):
+    opm_trace = opm_tg.getZTrace( i , j )
+    ecl_trace = ecl_tg.getZTrace( i , j )
+
+    print "OPM: %s" % opm_trace
+    print "ECL: %s" % ecl_trace
+    print
+    
+
+#-----------------------------------------------------------------
+
+if len(sys.argv) < 3:
+    sys.exit("example.py  opm_trans.json  eclipse_trans.json")
+
+opm_tg = TransGraph.load(sys.argv[1])
+ecl_tg = TransGraph.load(sys.argv[2])
+
+
+direction_example( opm_tg , ecl_tg )
+cell_example( opm_tg , ecl_tg )
+count_example( opm_tg , ecl_tg )
+
+
+xtrace_example( opm_tg , ecl_tg , 20 ,20)
+ytrace_example( opm_tg , ecl_tg , 10 ,10)
+
+ztrace_example( opm_tg , ecl_tg , 10 ,10)
+ztrace_example( opm_tg , ecl_tg , 10 ,20)
+ztrace_example( opm_tg , ecl_tg , 20 ,10)
+ztrace_example( opm_tg , ecl_tg , 20 ,20)
+ztrace_example( opm_tg , ecl_tg , 30 ,70)
+ztrace_example( opm_tg , ecl_tg , 30 ,60)

examples/python/trans_graph.py

@@ -0,0 +1,198 @@
+# The opm_init_check cpp program will produce two json files which
+# should describe the transmissibility graph. The structure of the
+# main chunk of data is a list:
+#
+#
+# ((i1,j1,k1) , (((i2,j2,k2), T12) , ((i3,j3,k3) , T13))),    
+# ((iA,jA,kA) , (((iB,jB,kB), TAB) , ((iC,jC,kC) , TAC))),
+# ....
+#
+#
+# Cell(i1,j1,k1) is connected to cells (i2,j2,k2) and (i3,j3,k3)
+# respectively, with transmissibility T12 and T13
+# respectively. Furthermore cell (iA,iB,iC) is connected to cells
+# (iB,jB,kB) and (iC,jC,kC) with transmissibilty TAB and TAC
+# respectively.
+
+
+import json
+
+
+
+
+class Connection(object):
+    """
+    The connection class holds connection information for one cell;
+    including the i,j,k of the cell and the Transmissibility of connection.
+    """
+
+    def __init__(self , i1, j1 , k1 , i2 , j2 , k2 , T):
+        self.i = i2 
+        self.j = j2
+        self.k = k2
+        self.T = T
+
+        dx = abs(i1 - i2)
+        dy = abs(j1 - j2)
+        dz = abs(k1 - k2)
+
+        if dx == 1 and dy == 0 and dz == 0:
+            self.dir = "X"
+        elif dx == 0 and dy == 1 and dz == 0:
+            self.dir = "Y"
+        elif dx == 0 and dy == 0 and dz == 1:
+            self.dir = "Z"
+        else:
+            self.dir = "NNC"
+
+
+    def __str__(self):
+        return "<%d,%d,%d>(T:%g)" % (self.i , self.j , self.k , self.T)
+
+
+
+class CellConnections(object):
+
+    def __init__(self , i,j,k):
+        self.i = i
+        self.j = j
+        self.k = k
+        self.connection_list = []
+        self.connection_map = {}
+        
+
+    def __getitem__(self , index):
+        if isinstance(index,int):
+            return self.connection_list[index]
+        else:
+            return self.connection_map[index]
+
+    def __len__(self):
+        return len(self.connection_list)
+
+    def has_key(self , dir_key):
+        return self.connection_map.has_key(dir_key)
+
+
+
+    def addConnections(self , connections):
+        for ijk,T in connections:
+            new_conn = Connection( self.i , self.j , self.k , ijk[0] , ijk[1] , ijk[2] , T)
+            self.connection_list.append( new_conn )
+            if new_conn.dir == "NNC":
+                if not self.connection_map.has_key("NNC"):
+                    self.connection_map["NNC"] = []
+                self.connection_map["NNC"].append( new_conn )
+            else:
+                self.connection_map[new_conn.dir] = new_conn
+
+
+    
+    def __nonzero__(self):
+        if len(self.connection_list) > 0:
+            return True
+        else:
+            return False
+
+
+    def connectsWith(self, i , j , k):
+        for conn in self.connection_list:
+            if conn.i == i and conn.j == j and conn.k == k:
+                return True
+            else:
+                return False
+
+
+    def __str__(self):
+        return "<%d,%d,%d> : %s" % (self.i , self.j , self.k , [ conn.__str__() for conn in self.connection_list ])
+
+
+
+class TransGraph(object):
+
+    def __init__(self , nx , ny , nz):
+        self.nx = nx
+        self.ny = ny
+        self.nz = nz
+    
+        self.cell_connections = [ None ] * nx * ny * nz
+
+
+    def __getitem__(self, index):
+        if isinstance(index , tuple):
+            g = index[0] + index[1] * self.nx + index[2]*self.nx*self.ny
+        else:
+            g = index
+
+        connections = self.cell_connections[g]
+        if connections is None:
+            k = g / (self.nx * self.ny)
+            j = (g - k * self.nx * self.ny) / self.nx
+            i = g - k * self.nx * self.ny - j * self.nx
+            self.cell_connections[g] = CellConnections( i,j,k )
+        
+        return self.cell_connections[g]
+        
+
+
+
+    def addCell(self , ijk , new_connections):
+        g = ijk[0] + ijk[1] * self.nx + ijk[2]*self.nx*self.ny
+        connections = self[g]
+        connections.addConnections( new_connections )
+
+
+
+    def getZTrace(self , i , j):
+        trace = []
+        for k in range(self.nz):
+            cell = self[i,j,k]
+            if cell.has_key("Z"):
+                trace.append( cell["Z"].T )
+            else:
+                trace.append( 0 )
+                
+        return trace
+
+
+    def getXTrace(self , j , k):
+        trace = []
+        for i in range(self.nx):
+            cell = self[i,j,k]
+            if cell.has_key("X"):
+                trace.append( cell["X"].T )
+            else:
+                trace.append( 0 )
+                
+        return trace
+
+
+    def getYTrace(self , i , k):
+        trace = []
+        for j in range(self.ny):
+            cell = self[i,j,k]
+            if cell.has_key("Y"):
+                trace.append( cell["Y"].T )
+            else:
+                trace.append( 0 )
+
+        return trace
+
+
+        
+    
+    @staticmethod
+    def load(json_file):
+        with open(json_file) as fileH:
+            data = json.load( fileH )
+        
+        dims = data["dims"]
+        graph = data["graph"]
+        
+        trans_graph = TransGraph( dims[0] , dims[1] , dims[2])
+
+        for cell,connections in graph:
+            trans_graph.addCell( cell , connections )
+        
+        return trans_graph
+