// ===========================================================================
//
//       Filename:  EMSchur3DBase.cpp
//
//        Created:  09/20/2013 04:53:40 PM
//       Compiler:  Tested with g++, icpc, and MSVC 2010
//
//         Author:  Trevor Irons (ti)
//
//   Organisation:  University of Utah
//
//          Email:  Trevor.Irons@utah.edu
//
// ===========================================================================

/**
  @file
  @author   Trevor Irons
  @date     09/20/2013
 **/

#include "EMSchur3DBase.h"

typedef Eigen::Triplet<Lemma::Complex> Tc;
typedef Eigen::Triplet<Lemma::Real> Tr;

// Moved to header
//#define UPPER 0  // LOWER WAS 0
//#define LOWER 1  // 1=true, 0=false

//#define FINITEVOLUME 1
#define FINITEDIFFERENCE 1

namespace Lemma {


    // ====================  FRIEND METHODS  =====================

    std::ostream &operator<<(std::ostream &stream, const EMSchur3DBase &ob) {
        stream << ob.Serialize()  << "\n";
        return stream;
    }


    // ====================  LIFECYCLE     =======================

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  EMSchur3DBase
    // Description:  constructor (protected)
    //--------------------------------------------------------------------------------------
//     std::shared_ptr<EMSchur3DBase> EMSchur3DBase::NewSP() {
//         return std::make_shared<EMSchur3DBase>( ctor_key() );
// 	}

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  EMSchur3DBase
    // Description:  constructor (protected)
    //--------------------------------------------------------------------------------------
    EMSchur3DBase::EMSchur3DBase ( const ctor_key& key ) : LemmaObject( key ),
        Grid(nullptr),
        //Survey(nullptr),
        LayModel(nullptr), Cvec(nullptr),
        ResFile("source"), sigma(nullptr), sigmap(nullptr)
    {
    }  // -----  end of method EMSchur3DBase::EMSchur3DBase  (constructor)  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  EMSchur3DBase
    // Description:  constructor (protected)
    //--------------------------------------------------------------------------------------
    EMSchur3DBase::EMSchur3DBase ( const YAML::Node& node, const ctor_key& key ) : LemmaObject( key ),
        Grid(nullptr),
        //Survey(nullptr),
        LayModel(nullptr), Cvec(nullptr),
        ResFile("source"), sigma(nullptr), sigmap(nullptr)
    {
    }  // -----  end of method EMSchur3DBase::~EMSchur3DBase  (destructor)  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  ~EMSchur3DBase
    // Description:  destructor (protected)
    //--------------------------------------------------------------------------------------
    EMSchur3DBase::~EMSchur3DBase ( ) {
        if (sigma) Delete3DScalar(sigma);
        if (sigmap) Delete3DScalar(sigmap);
        if (Cvec) delete [] Cvec;
        //if (CvecRe) delete [] CvecRe;
    }  // -----  end of method EMSchur3DBase::~EMSchur3DBase  (destructor)  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3D
    //      Method:  Serialize
    //--------------------------------------------------------------------------------------
    YAML::Node EMSchur3DBase::Serialize (  ) const {
        YAML::Node node = LemmaObject::Serialize();
        node.SetTag( this->GetName() );
        node["EMSchur3D_VERSION"] = EMSCHUR3D_VERSION;
        return node;
    }		// -----  end of method EMSchur3D::Serialize  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  BuildC
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::BuildC ( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw) {

        Cvec[iw].resize( unx+uny+unz , unx+uny+unz );

#if LOWER && UPPER
        Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 7));   // Whole
#else
        Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 4)); // Upper/Lower
#endif

        //Cvec_s.resize( idx. )
        //CMMvec[iw].resize( unx+uny+unz , unx+uny+unz );
        //CMMvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 1));

        Real omega = Omegas[iw];
        std::cout << "Building C_" << iw << std::endl;

        //Complex hsig(0);
        Real hsig(0);
        Real EPS(1e-24);
        Real EPS2(1e-18);

        // book keeping
        int ir = 0;
        int ic = 0;

        // LAPL{Ax} + i omega mu sigmax
        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny;   ++iy) {
        for (int ix=0; ix<nx+1; ++ix) {

            // Calculate 1/2 step values on staggered grid
            Real hzlo(0);
            Real hzhi(0);
            if (iz == 0) {
                hzlo = .5*hz[iz+1] + .5*hz[iz];                   // assume same grid spacing for boundary
                hzhi = .5*hz[iz+1] + .5*hz[iz];
            } else if (iz == nz-1) {
                hzlo = .5*hz[iz-1] + .5*hz[iz];
                hzhi = .5*hz[iz-1] + .5*hz[iz];                   // assume same grid spacing for boundary
            } else {
                hzlo = .5*hz[iz-1] + .5*hz[iz  ];
                hzhi = .5*hz[iz  ] + .5*hz[iz+1];
            }
            Real scz = 2./(hzlo+hzhi);

            // Calculate 1/2 step values on staggered grid
            Real hylo(0);
            Real hyhi(0);
            if (iy == 0) {
                hylo = .5*hy[iy+1] + .5*hy[iy];                   // assume same grid spacing for boundary
                hyhi = .5*hy[iy+1] + .5*hy[iy];
            } else if (iy == ny-1) {
                hylo = .5*hy[iy-1] + .5*hy[iy];
                hyhi = .5*hy[iy-1] + .5*hy[iy];                   // assume same grid spacing for boundary
            } else {
                hylo = .5*hy[iy-1] + .5*hy[iy  ];
                hyhi = .5*hy[iy  ] + .5*hy[iy+1];
            }
            Real scy = 2./(hylo+hyhi);

            /////////////////////////////////////////////////////////////////////////
            // Diagonal entry
            // Harmonically average sigmax
            if (ix == 0) {
                hsig = sigma[ix][iy][iz];
            } else if (ix == nx) {
                hsig = sigma[ix-1][iy][iz];
            } else {
                hsig = ((hx[ix]+hx[ix-1])/2.) *
                        (1. / ( (hx[ix  ] / (2.*sigma[ix  ][iy][iz] + EPS))  +
                              (  hx[ix-1] / (2.*sigma[ix-1][iy][iz] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            //hsig += Complex(0, omega*EPSILON0);

            Real Sum(0);
            Real scx(0);
            if (ix == 0) { // x- bdry
                scx = 2./(hx[ix]+hx[ix]);
                if (DirichletXLO) {
                    /* Dirichlet bdry condition on A */
                    Sum = scz/hzlo + scz/hzhi + scy/hylo + scy/hyhi + 2.*scx/hx[ix];
                } else {
                    /* Neumann bdry */
                    Sum = scz/hzlo + scz/hzhi + scy/hylo + scy/hyhi + scx/hx[ix];
                }
            } else if (ix == nx) {  // x+ bdry
                scx = 2./(hx[ix-1]+hx[ix-1]);
                if (DirichletXHI) {
                    /* Dirichlet bdry condition on A */
                    Sum = scz/hzlo + scz/hzhi + scy/hylo + scy/hyhi + 2.*scx/hx[ix-1];
                } else {
                    /* Neumann bdry */
                    Sum = scz/hzlo + scz/hzhi + scy/hylo + scy/hyhi + scx/hx[ix-1];
                }
            } else {
                scx = 2./(hx[ix]+hx[ix-1]);
                Sum = scz/hzlo + scz/hzhi + scy/hylo + scy/hyhi + scx/hx[ix] + scx/hx[ix-1];
            }

            /////////////////////////////////////////
            // minus side off diagonal terms
#if LOWER
            // Third Off Diagonal
            if (iz!=0)         Cvec[iw].insert(ir, ic-ny*(nx+1)) = Complex(-scz/hzlo, 0);

            // Second Off Diagonal
            if (iy!=0)         Cvec[iw].insert(ir, ic-(nx+1)) = Complex(-scy/hylo, 0);

            // First  Off Diagonal
            if (ix!=0)         Cvec[iw].insert(ir, ic-1) = Complex(-scx/hx[ix-1], 0);
#endif
            ////////////////////////////////////////////
            // Diagonal Term
            Cvec[iw].insert(ir,ic) = Complex(Sum, 0) + Complex(0,1)*omega*MU0*hsig;
            //Cvec[iw].insert(ir,ic) =        Complex(Sum, omega*MU0*hsig);
            //CMMvec[iw].insert(ir,ic) =  1./ Complex(Sum, omega*MU0*hsig);

            ////////////////////////////////////////////////////////////////////////
            // plus side off diagonal terms
#if UPPER
            // First Off Diagonal
            if (ix!=nx)         Cvec[iw].insert(ir, ic+1) = Complex(-scx/hx[ix], 0);

            // Second Off Diagonal
            if (iy!=ny-1)       Cvec[iw].insert(ir, ic+(nx+1)) = Complex(-scy/hyhi, 0);

            // Third Off Diagonal
            if (iz!=nz-1)       Cvec[iw].insert(ir, ic+ny*(nx+1)) = Complex(-scz/hzhi, 0);
#endif
            ++ir;
            ++ic;
        }
        }
        }
        assert(ic == unx);

        // LAPL{Ay} + i omega mu sigmay
        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny+1; ++iy) {
        for (int ix=0; ix<nx;   ++ix) {

            // Calculate 1/2 step values on staggered grid
            Real hxlo(0);
            Real hxhi(0);
            if (ix == 0) {
                hxlo = .5*hx[ix+1] + .5*hx[ix];                   // assume same grid spacing for boundary
                hxhi = .5*hx[ix+1] + .5*hx[ix];
            } else if (ix == nx-1) {
                hxlo = .5*hx[ix-1] + .5*hx[ix];
                hxhi = .5*hx[ix-1] + .5*hx[ix];                   // assume same grid spacing for boundary
            } else {
                hxlo = .5*hx[ix-1] + .5*hx[ix  ];
                hxhi = .5*hx[ix  ] + .5*hx[ix+1];
            }
            Real scx = 2./(hxlo+hxhi);

            // Calculate 1/2 step values on staggered grid
            Real hzlo(0);
            Real hzhi(0);
            if (iz == 0) {
                hzlo = .5*hz[iz+1] + .5*hz[iz];                   // assume same grid spacing for boundary
                hzhi = .5*hz[iz+1] + .5*hz[iz];
            } else if (iz == nz-1) {
                hzlo = .5*hz[iz-1] + .5*hz[iz];
                hzhi = .5*hz[iz-1] + .5*hz[iz];                   // assume same grid spacing for boundary
            } else {
                hzlo = .5*hz[iz-1] + .5*hz[iz  ];
                hzhi = .5*hz[iz  ] + .5*hz[iz+1];
            }
            Real scz = 2./(hzlo+hzhi);

            // Harmonically average sigmay
            if (iy == 0) {
                hsig = sigma[ix][iy][iz];
            } else if (iy == ny) {
                hsig = sigma[ix][iy-1][iz];
            } else {
                hsig = ((hy[iy]+hy[iy-1])/2.) *
                        (1. / ( (hy[iy  ] / (2.*sigma[ix][iy  ][iz] + EPS))  +
                              (  hy[iy-1] / (2.*sigma[ix][iy-1][iz] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            //hsig += Complex(0, omega*EPSILON0);

            Real Sum(0);
            Real scy(0);
            if (iy == 0) {  // y- bdry
                scy = 2./(hy[iy]+hy[iy]);
                if (DirichletYLO) {
                    /* Dirichlet bdry condition on A */
                    Sum = scz/hzlo + scz/hzhi + scx/hxlo + scx/hxhi + 2.*scy/hy[iy];
                } else {
                    /* Neumann bdry */
                    Sum = scz/hzlo + scz/hzhi + scx/hxlo + scx/hxhi + scy/hy[iy];
                }
            } else if (iy == ny) {  // y+ bdry
                scy = 2./(hy[iy-1]+hy[iy-1]);
                if (DirichletYHI) {
                    /* Dirichlet bdry condition on A */
                    Sum = scz/hzlo + scz/hzhi + scx/hxlo + scx/hxhi + 2.*scy/hy[iy-1];
                } else {
                    /* Neumann bdry */
                    Sum = scz/hzlo + scz/hzhi + scx/hxlo + scx/hxhi + scy/hy[iy-1];
                }
            } else {
                scy = 2./(hy[iy]+hy[iy-1]);
                Sum = scz/hzlo + scz/hzhi + scx/hxlo + scx/hxhi + scy/hy[iy-1] + scy/hy[iy];
            }
#if LOWER
            // Third Off Diagonal
            if (iz!=0)              Cvec[iw].insert(ir,ic-(ny+1)*(nx)) = Complex(-scz/hzlo, 0);

            // Second Off Diagonal
            if (iy!=0)              Cvec[iw].insert(ir,ic-(nx)) = Complex(-scy/hy[iy-1], 0);

            // First Off Diagonal
            if (ix!=0)              Cvec[iw].insert(ir,ic-1) = Complex(-scx/hxlo, 0);
#endif
            // Diagonal Term
            Cvec[iw].insert(ir,ic) =       Complex(Sum, 0) + Complex(0,1)*omega*MU0*hsig;
            //Cvec[iw].insert(ir,ic) =        Complex(Sum, omega*MU0*hsig);
            //CMMvec[iw].insert(ir,ic) =  1./ Complex(Sum, omega*MU0*hsig);
#if UPPER
            // First Off Diagonal
            if (ix!=nx-1)           Cvec[iw].insert(ir,ic+1) = Complex(-scx/hxhi, 0);

            // Second Off Diagonal
            if (iy!=ny)             Cvec[iw].insert(ir,ic+(nx)) = Complex(-scy/hy[iy], 0);

            // Third Off Diagonal
            if (iz!=nz-1)           Cvec[iw].insert(ir,ic+(ny+1)*(nx)) = Complex(-scz/hzhi, 0);
#endif
            ++ir;
            ++ic;
        }
        }
        }
        assert(ic == unx+uny);

        // LAPL{Az} + i omega mu sigmaz
        for (int iz=0; iz<nz+1; ++iz) {
        for (int iy=0; iy<ny;   ++iy) {
        for (int ix=0; ix<nx;   ++ix) {

            // Calculate 1/2 step values on staggered grid
            Real hylo(0);
            Real hyhi(0);
            if (iy == 0) {
                hylo = .5*hy[iy+1] + .5*hy[iy];                   // assume same grid spacing for boundary
                hyhi = .5*hy[iy+1] + .5*hy[iy];
            } else if (iy == ny-1) {
                hylo = .5*hy[iy-1] + .5*hy[iy];
                hyhi = .5*hy[iy-1] + .5*hy[iy];                   // assume same grid spacing for boundary
            } else {
                hylo = .5*hy[iy-1] + .5*hy[iy  ];
                hyhi = .5*hy[iy  ] + .5*hy[iy+1];
            }
            Real scy = 2./(hylo+hyhi);

            // Calculate 1/2 step values on staggered grid
            Real hxlo(0);
            Real hxhi(0);
            if (ix == 0) {
                hxlo = .5*hx[ix+1] + .5*hx[ix];                   // assume same grid spacing for boundary
                hxhi = .5*hx[ix+1] + .5*hx[ix];
            } else if (ix == nx-1) {
                hxlo = .5*hx[ix-1] + .5*hx[ix];
                hxhi = .5*hx[ix-1] + .5*hx[ix];                   // assume same grid spacing for boundary
            } else {
                hxlo = .5*hx[ix-1] + .5*hx[ix  ];
                hxhi = .5*hx[ix  ] + .5*hx[ix+1];
            }
            Real scx = 2./(hxlo+hxhi);

            // Harmonically average sigmaz
            if (iz==0) {
                hsig = sigma[ix][iy][iz];
            } else if (iz == nz) {
                hsig = sigma[ix][iy][iz-1];
            } else {
                hsig = ((hz[iz]+hz[iz-1])/2.) *
                        (1. / ( (hz[iz  ] / (2.*sigma[ix][iy][iz  ] + EPS))  +
                              (  hz[iz-1] / (2.*sigma[ix][iy][iz-1] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            //hsig += Complex(0, omega*EPSILON0);

            Real Sum(0);
            Real scz(0);
            if (iz == 0) { // z- bdry
                scz = 2./(hz[iz]+hz[iz]);
                if (DirichletZLO) {
                    /* Dirichlet bdry condition on A */
                    Sum = scx/hxlo + scx/hxhi + scy/hylo + scy/hyhi + 2.*scz/hz[iz];
                } else {
                    // Neumann bdry
                    Sum = scx/hxlo + scx/hxhi + scy/hylo + scy/hyhi + scz/hz[iz];
                }
            } else if (iz == nz) { // z+ bdry
                scz = 2./(hz[iz-1]+hz[iz-1]);
                if (DirichletZHI) {
                    // Dirichlet bdry condition on A
                    Sum = scx/hxlo + scx/hxhi + scy/hylo + scy/hyhi + 2.*scz/hz[iz-1];
                } else {
                    // Neumann bdry
                    //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz-1];
                    Sum = scx/hxlo + scx/hxhi + scy/hylo + scy/hyhi + scz/hz[iz-1];
                }
            } else {
                scz = 2./(hz[iz]+hz[iz-1]);
                Sum = scx/hxlo + scx/hxhi + scy/hylo + scy/hyhi + scz/hz[iz-1] + scz/hz[iz];
            }
#if LOWER
            // Third Off Diagonal
            if (iz!=0)             Cvec[iw].insert(ir,ic-ny*(nx)) = Complex(-scz/hz[iz-1], 0);

            // Second Off Diagonal
            if (iy!=0)           Cvec[iw].insert(ir,ic-(nx)) = Complex(-scy/hylo, 0);

            // First Off Diagonal
            if (ix!=0)           Cvec[iw].insert(ir,ic-1) = Complex(-scx/hxlo, 0);
#endif
            // Diagonal Term
            Cvec[iw].insert(ir,ic) =       Complex(Sum, 0) + Complex(0,1)*omega*MU0*hsig;
            //Cvec[iw].insert(ir,ic) =         Complex(Sum, omega*MU0*hsig);
            //CMMvec[iw].insert(ir,ic) =  1. / Complex(Sum, omega*MU0*hsig);
#if UPPER
            // First Off Diagonal
            if (ix!=nx-1)           Cvec[iw].insert(ir,ic+1) = Complex(-scx/hx[ix], 0);

            // Second Off Diagonal
            if (iy!=ny-1)           Cvec[iw].insert(ir,ic+(nx)) = Complex(-scy/hy[iy], 0);

            // Third Off Diagonal
            if (iz!=nz)             Cvec[iw].insert(ir,ic+ny*(nx)) = Complex(-scz/hz[iz], 0);
#endif
            ++ir;
            ++ic;
        }
        }
        }
        assert(ic == unx+uny+unz);
        //Cvec[iw].makeCompressed();
        Cvec[iw].makeCompressed();
        //CMMvec[iw].makeCompressed();
        //CsymVec[iw] = Cvec[iw].selfadjointView<Eigen::Upper>();
        return ;
    }		// -----  end of method EMSchur3DBase::BuildC  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  BuildCReal
    //--------------------------------------------------------------------------------------
    /*
    void EMSchur3DBase::BuildCReal ( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw) {

        int CI = unx+uny+unz;
        CvecRe[iw].resize( 2*CI, 2*CI );       // twice the size, but Real

        std::vector<Tr> triplets;
        triplets.reserve( 9*CI );

        Real omega = Omegas[iw];
        std::cout << "Building real C_" << iw << std::endl;

        // LAPL{Ax} + i omega mu sigmax
        int ir = 0;
        int ic = 0;
        Real hsig(0);
        Real EPS(1e-24);
        Real EPS2(1e-18);

        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny;   ++iy) {
        for (int ix=0; ix<nx+1; ++ix) {

            // Calculate 1/2 step values on staggered grid
            Real alo_ihz2(0);
            Real ahi_ihz2(0);
            if (iz == 0) {
                ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
                alo_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
            } else if (iz == nz-1) {
                alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
                ahi_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
            } else {
                alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
                ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
            }

            Real alo_ihy2(0);  // half step low jump in y
            Real ahi_ihy2(0);  // half step high jump in y
            if (iy == 0) {
                ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
                alo_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
            } else if (iy == ny-1) {
                alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
                ahi_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
            } else {
                alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
                ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
            }

            /////////////////////////////////////////////////////////////////////////
            // Diagonal entry
            // Harmonically average sigmax
            if (ix == 0) {
                hsig = sigma[ix][iy][iz];
            } else if (ix == nx) {
                hsig = sigma[ix-1][iy][iz];
            } else {
                hsig = ((hx[ix]+hx[ix-1])/2.) *
                        (1. / ( (hx[ix  ] / (2.*sigma[ix  ][iy][iz] + EPS))  +
                              (  hx[ix-1] / (2.*sigma[ix-1][iy][iz] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            //hsig += Complex(0, omega*EPSILON0);

            Real Sum(0);
            if (ix == 0) {
                // Dirichlet bdry condition on A
                // Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + 2.*ihx2[ix  ];
                // Neumann bdry
                Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix  ];
            } else if (ix == nx) {
                // Dirichlet bdry condition on A
                // Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + 2.*ihx2[ix-1];
                // Neumann bdry
                Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix-1];
            } else {
                Sum = alo_ihy2+ahi_ihy2 + alo_ihz2+ahi_ihz2 + ihx2[ix] + ihx2[ix-1];
            }

            ////////////////////////////////////////////
            // Diagonal Term
            triplets.push_back( Tr(ir,    ic   ,  Sum) );
            triplets.push_back( Tr(ir+CI, ic+CI, -Sum) );

            ////////////////////////////////////////////////////////////////////////
            // plus side off diagonal terms
            // First Off Diagonal
            if (ix!=nx) {
                triplets.push_back( Tr(ir,    ic+1   , -ihx2[ix]) );
                triplets.push_back( Tr(ir+CI, ic+1+CI,  ihx2[ix]) );
            }
            // Second Off Diagonal
            if (iy!=ny-1) {
                triplets.push_back( Tr(ir   , ic+(nx+1)   , -ahi_ihy2) );
                triplets.push_back( Tr(ir+CI, ic+(nx+1)+CI,  ahi_ihy2) );
            }
            // Third Off Diagonal
            if (iz!=nz-1) {
                triplets.push_back( Tr(ir   , ic+ny*(nx+1)   , -ahi_ihz2) );
                triplets.push_back( Tr(ir+CI, ic+ny*(nx+1)+CI,  ahi_ihz2) );
            }

            ////////////////////////////////////////////////////////////////////////
            // imaginary part
            triplets.push_back( Tr(ir, ic+CI, omega*MU0*hsig) );

            ++ir;
            ++ic;
        }
        }
        }
        assert(ic == unx);

        // LAPL{Ay} + i omega mu sigmay
        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny+1; ++iy) {
        for (int ix=0; ix<nx;   ++ix) {

            // Calculate 1/2 step values on staggered grid
            Real alo_ihz2(0);
            Real ahi_ihz2(0);
            if (iz == 0) {
                ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
                alo_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
            } else if (iz == nz-1) {
                alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
                ahi_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
            } else {
                alo_ihz2 = std::pow(.5*hz[iz-1] + .5*hz[iz], -2);
                ahi_ihz2 = std::pow(.5*hz[iz+1] + .5*hz[iz], -2);
            }

            // Calculate 1/2 step values on staggered grid
            Real alo_ihx2(0);
            Real ahi_ihx2(0);
            if (ix == 0) {
                ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
                alo_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
            } else if (ix == nx-1) {
                alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
                ahi_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
            } else {
                alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
                ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
            }

            // Harmonically average sigmay
            if (iy == 0) {
                hsig = sigma[ix][iy][iz];
            } else if (iy == ny) {
                hsig = sigma[ix][iy-1][iz];
            } else {
                hsig = ((hy[iy]+hy[iy-1])/2.) *
                        (1. / ( (hy[iy  ] / (2.*sigma[ix][iy  ][iz] + EPS))  +
                              (  hy[iy-1] / (2.*sigma[ix][iy-1][iz] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            //hsig += Complex(0, omega*EPSILON0);

            Real Sum(0);
            if (iy == 0) {
                // Dirichlet bdry condition on A
                //Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + 2.*ihy2[iy  ];
                // Neumann bdry
                Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy  ];
            } else if (iy == ny) {
                // Dirichlet bdry condition on A
                //Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + 2.*ihy2[iy-1];
                // Neumann bdry
                Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy-1];
            } else {
                Sum = alo_ihx2+ahi_ihx2 + alo_ihz2+ahi_ihz2 + ihy2[iy] + ihy2[iy-1];
            }

            //////////////////////////////////////////////////////////////////////////////
            // Diagonal Term
            triplets.push_back( Tr(ir,    ic   ,  Sum) );
            triplets.push_back( Tr(ir+CI, ic+CI, -Sum) );

            // First Off Diagonal
            if (ix!=nx-1) {
                triplets.push_back( Tr(ir,    ic+1   , -ahi_ihx2));
                triplets.push_back( Tr(ir+CI, ic+1+CI,  ahi_ihx2));
            }

            // Second Off Diagonal
            if (iy!=ny) {
                triplets.push_back( Tr(ir   , ic+(nx)   , -ihy2[iy]) );
                triplets.push_back( Tr(ir+CI, ic+(nx)+CI,  ihy2[iy]) );
            }

            // Third Off Diagonal
            if (iz!=nz-1) {
                triplets.push_back( Tr(ir   , ic+(ny+1)*(nx)   , -ahi_ihz2));
                triplets.push_back( Tr(ir+CI, ic+(ny+1)*(nx)+CI,  ahi_ihz2));
            }

            ///////////////////////////////////////////////////////////////////////
            // imaginary part
            triplets.push_back( Tr(ir, ic+CI, omega*MU0*hsig) );

            ++ir;
            ++ic;
        }
        }
        }
        assert(ic == unx+uny);

        // LAPL{Az} + i omega mu sigmaz
        for (int iz=0; iz<nz+1; ++iz) {
        for (int iy=0; iy<ny;   ++iy) {
        for (int ix=0; ix<nx;   ++ix) {

            // Calculate 1/2 step values on staggered grid
            Real alo_ihx2(0);
            Real ahi_ihx2(0);
            if (ix == 0) {
                ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
                alo_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
            } else if (ix == nx-1) {
                alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
                ahi_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
            } else {
                alo_ihx2 = std::pow(.5*hx[ix-1] + .5*hx[ix], -2);
                ahi_ihx2 = std::pow(.5*hx[ix+1] + .5*hx[ix], -2);
            }

            // Calculate 1/2 step values on staggered grid
            Real alo_ihy2(0);
            Real ahi_ihy2(0);
            if (iy == 0) {
                ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
                alo_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
            } else if (iy == ny-1) {
                alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
                ahi_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
            } else {
                alo_ihy2 = std::pow(.5*hy[iy-1] + .5*hy[iy], -2);
                ahi_ihy2 = std::pow(.5*hy[iy+1] + .5*hy[iy], -2);
            }

            // Harmonically average sigmaz
            if (iz==0) {
                hsig = sigma[ix][iy][iz];
            } else if (iz == nz) {
                hsig = sigma[ix][iy][iz-1];
            } else {
                hsig = ((hz[iz]+hz[iz-1])/2.) *
                        (1. / ( (hz[iz  ] / (2.*sigma[ix][iy][iz  ] + EPS))  +
                              (  hz[iz-1] / (2.*sigma[ix][iy][iz-1] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            //hsig += Complex(0, omega*EPSILON0);

            Real Sum(0);
            if (iz == 0) {
                // Dirichlet bdry condition on A
                //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + 2.*ihz2[iz  ];
                // Neumann bdry
                Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz  ];
            } else if (iz == nz) {
                // Dirichlet bdry condition on A
                //Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + 2.*ihz2[iz-1];
                // Neumann bdry
                Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz-1];
            } else {
                Sum = alo_ihx2 + ahi_ihx2 + alo_ihy2 + ahi_ihy2 + ihz2[iz] + ihz2[iz-1];
            }

            //////////////////////////////////////////////////////////////////
            // Diagonal Term
            triplets.push_back( Tr(ir   , ic   ,  Sum) );
            triplets.push_back( Tr(ir+CI, ic+CI, -Sum) );

            // First Off Diagonal
            if (ix!=nx-1) {
                triplets.push_back( Tr(ir   , ic+1   , -ahi_ihx2) );
                triplets.push_back( Tr(ir+CI, ic+1+CI,  ahi_ihx2) );
            }
            // Second Off Diagonal
            if (iy!=ny-1) {
                triplets.push_back( Tr(ir   , ic+(nx)   , -ahi_ihy2) );
                triplets.push_back( Tr(ir+CI, ic+(nx)+CI,  ahi_ihy2) );
            }
            // Third Off Diagonal
            if (iz!=nz) {
                triplets.push_back( Tr(ir   , ic+ny*(nx)   , -ihz2[iz]) );
                triplets.push_back( Tr(ir+CI, ic+ny*(nx)+CI,  ihz2[iz]) );
            }

            //////////////////////////////////////////////////////
            // imaginary part
            triplets.push_back( Tr(ir, ic+CI, omega*MU0*hsig) );

            ++ir;
            ++ic;
        }
        }
        }
        assert(ic == unx+uny+unz);
        CvecRe[iw].setFromTriplets(triplets.begin(), triplets.end()) ;
        CvecRe[iw].makeCompressed();

        return ;
    }		// -----  end of method EMSchur3DBase::BuildCReal  -----
    */

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  SetResFileName
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::SetResFileName ( const std::string& fname ) {
        ResFile = fname;
        return ;
    }		// -----  end of method EMSchur3DBase::SetResFileName  -----


    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  BuildCPreconditioner
    //--------------------------------------------------------------------------------------
//     void EMSchur3DBase::BuildCPreconditioner ( const int& iw ) {
//         std::cout << "Building preconditioner for C_" << iw << std::endl;
//         CMvec[iw].setDroptol(1e-4); // 1e-4 previously (tested value)
//         CMvec[iw].setFillfactor(10);
//         Eigen::SparseMatrix<Complex>  Csym;
//         Csym = Cvec[iw].selfadjointView<Eigen::Upper>();
//         CMvec[iw].compute( Csym );
//         return ;
//     }		// -----  end of method EMSchur3DBase::BuildCPreconditioner  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  BuildD
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::BuildD ( ) {
        D.resize(uns, unx+uny+unz );
        //D.reserve(Eigen::VectorXi::Constant(unx+uny+unz, 6));

        std::vector<Tc> triplets;
        triplets.reserve( 2*unx + 2*uny + 2*unz );

        std::cout << "Building D...";
        int ic = 0;
        int ir = 0;

        ////////////////////////////////
        // Ax dx
        for (int iz=0; iz<nz; ++iz) {
        for (int iy=0; iy<ny; ++iy) {
            ++ic; // TRICKY
            for (int ix=0; ix<nx; ++ix) {
                //D.insert(ir, ic-1) = Complex(-ihx[ix], 0);
                //D.insert(ir, ic  ) = Complex( ihx[ix], 0);
#ifdef FINITEDIFFERENCE
                //triplets.push_back( Tc(ir, ic-1, Complex(-sc*ihx[ix], 0)) );
                //triplets.push_back( Tc(ir, ic  , Complex( sc*ihx[ix], 0)) );
                triplets.push_back( Tc(ir, ic-1, Complex(-1.0/hx[ix], 0)) );
                triplets.push_back( Tc(ir, ic  , Complex( 1.0/hx[ix], 0)) );
#elif FINITEVOLUME
                triplets.push_back( Tc(ir, ic-1, Complex(-1, 0)) );
                triplets.push_back( Tc(ir, ic  , Complex( 1, 0)) );
#endif
                ++ir;
                ++ic;
            }
        }
        }
        assert (ic==unx);

        ///////////////////////////////
        // Ay dy
        ir = 0;
        for (int iz=0; iz<nz; ++iz) {
            for (int iy=0; iy<ny; ++iy) {
            for (int ix=0; ix<nx; ++ix) {
                //D.insert(ir, ic   ) = Complex(-ihy[iy], 0);
                //D.insert(ir, ic+nx) = Complex( ihy[iy], 0);
#ifdef FINITEDIFFERENCE
                //triplets.push_back( Tc(ir, ic   , Complex(-sc*ihy[iy], 0)) );
                //triplets.push_back( Tc(ir, ic+nx, Complex( sc*ihy[iy], 0)) );
                triplets.push_back( Tc(ir, ic   , Complex(-1.0/hy[iy], 0)) );
                triplets.push_back( Tc(ir, ic+nx, Complex( 1.0/hy[iy], 0)) );
#elif FINITEVOLUME
                triplets.push_back( Tc(ir, ic   , Complex(-1, 0)) );
                triplets.push_back( Tc(ir, ic+nx, Complex( 1, 0)) );
#endif
                ++ir;
                ++ic;
            }
            }
            ic += nx;
        }
        assert (ic==unx+uny);

        ///////////////////////////////
        // Az dz
        ir = 0;
        for (int iz=0; iz<nz; ++iz) {
        for (int iy=0; iy<ny; ++iy) {
        for (int ix=0; ix<nx; ++ix) {
            //D.insert(ir, ic      ) = Complex(-ihz[iz], 0);
            //D.insert(ir, ic+nx*ny) = Complex( ihz[iz], 0);
#ifdef FINITEDIFFERENCE
            //triplets.push_back( Tc(ir, ic      , Complex(-sc*ihz[iz], 0)) );
            //triplets.push_back( Tc(ir, ic+nx*ny, Complex( sc*ihz[iz], 0)) );
            triplets.push_back( Tc(ir, ic      , Complex(-1.0/hz[iz], 0)) );
            triplets.push_back( Tc(ir, ic+nx*ny, Complex( 1.0/hz[iz], 0)) );
#elif FINITEVOLUME
            triplets.push_back( Tc(ir, ic      , Complex(-1, 0)) );
            triplets.push_back( Tc(ir, ic+nx*ny, Complex( 1, 0)) );
#endif
            ++ir;
            ++ic;
        }
        }
        }
        assert (ic+nx*ny==unx+uny+unz);

        D.setFromTriplets(triplets.begin(), triplets.end());
        D.makeCompressed();
        std::cout << "Done!" << std::endl;
        return ;
    }		// -----  end of method EMSchur3DBase::BuildD  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  SetRectilinearGrid
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::SetRectilinearGrid ( std::shared_ptr<RectilinearGrid> GridPtr ) {

        Grid = GridPtr;

        // Do some convienience stuff, we call these so often it's nice to not have to use
        // the accessors every time. This is the only place these are set.
        nx = Grid->GetNx();
        ny = Grid->GetNy();
        nz = Grid->GetNz();

        hx = Grid->GetDx();
        hy = Grid->GetDy();
        hz = Grid->GetDz();

        unx  = (nx+1)*ny*nz;    // Number of Vectors x component
        uny  = nx*(ny+1)*nz;    // Number of Vectors y component
        unz  = nx*ny*(nz+1);    // Number of Vectors z component
        uns  = nx*ny*nz;        // On dual grid, number of scalars

        ihx = 1. / hx.array();
        ihy = 1. / hy.array();
        ihz = 1. / hz.array();

        ihx2 = 1. / (hx.array() * hx.array());
        ihy2 = 1. / (hy.array() * hy.array());
        ihz2 = 1. / (hz.array() * hz.array());

        return ;
    }		// -----  end of method EMSchur3DBase::SetRectilinearGrid  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  SetRectilinearGrid
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::SetRectilinearGrid ( vtkRectilinearGrid*  rGridPtr ) {

        // Convert to Lemma here


    }		// -----  end of method EMSchur3DBase::SetRectilinearGrid  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  Setup
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::Setup (  ) {

        ////////////////////////////////////
        // First set up grid stuff
        if (Cvec) delete [] Cvec;
        //if (CSolver) delete [] CSolver;
        Cvec    = new Eigen::SparseMatrix<Complex, Eigen::ColMajor> [Omegas.size()];

        //#ifdef LEMMAUSEOMP
        //#pragma omp parallel for schedule(static, 1)
        //#endif
        for (int iw=0; iw<Omegas.size(); ++iw) {
            std::cout << "Build C " << std::endl;
            BuildC(sigma, sigma, sigma, iw);
//             std::cout << "Build Re(C) " << std::endl;
//             BuildCReal(sigma, sigma, sigma, iw);
        }

        BuildCDirectSolver( ); // templated specialisation
        BuildD();

        return ;
    }		// -----  end of method EMSchur3DBase::Setup  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  SetLayeredEarthEM
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::SetLayeredEarthEM ( std::shared_ptr<LayeredEarthEM> LayModelin ) {
        LayModel = LayModelin;
        return ;
    }		// -----  end of method EMSchur3DBase::SetLayeredEarthEM  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  PrimaryField
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::PrimaryField ( std::shared_ptr<DipoleSource> source, std::shared_ptr<FieldPoints> dpoint ) {

        auto EmEarth = Lemma::EMEarth1D::NewSP();
            EmEarth->AttachDipoleSource(source);
            EmEarth->AttachLayeredEarthEM(LayModel);
            EmEarth->AttachFieldPoints(dpoint);
            EmEarth->SetFieldsToCalculate(Lemma::E);
            EmEarth->SetHankelTransformMethod(ANDERSON801);
            std::cout << "Primary Field Calculation" << std::endl;
            EmEarth->MakeCalc3();
        return ;
    }		// -----  end of method EMSchur3DBase::PrimaryField  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  FillPoints
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::FillPoints ( std::shared_ptr<FieldPoints> dpoint) {

        dpoint->SetNumberOfPoints(unx + uny + unz);
        int ip(0);
        Real tx, ty, tz;
        Real ox = Grid->GetOx();
        Real oy = Grid->GetOy();
        Real oz = Grid->GetOz();
        tz = oz;
        for (int iz=0; iz<nz;   ++iz) {
            ty = oy;
            for (int iy=0; iy<ny;   ++iy) {
                tx = ox - hx(0)/2.;
                for (int ix=0; ix<nx+1; ++ix) {
                    dpoint->SetLocation(ip, tx, ty, tz);
                    if (ix < nx) tx += hx(ix);
                    ++ip;
                }
                if (iy<ny-1) ty += .5*hy(iy) + .5*hy(iy+1);
            }
            if (iz<nz-1) tz += .5*hz(iz) + .5*hz(iz+1);
        }

        tz = oz;
        for (int iz=0; iz<nz;   ++iz) {
            ty = oy - hy(0)/2.;
            for (int iy=0; iy<ny+1; ++iy) {
                tx = ox;
                for (int ix=0; ix<nx;   ++ix) {
                    dpoint->SetLocation(ip, tx, ty, tz);
                    if (ix<nx-1) tx += .5*hx(ix) + .5*hx(ix+1);
                    ++ip;
                }
                if (iy < ny) ty += hy(iy);
            }
            if (iz<nz-1) tz += .5*hz(iz) + .5*hz(iz+1);
        }

        tz = oz - hz(0)/2.;
        for (int iz=0; iz<nz+1; ++iz) {
            ty = oy;
            for (int iy=0; iy<ny;   ++iy) {
                tx = ox;
                for (int ix=0; ix<nx;   ++ix) {
                    dpoint->SetLocation(ip, tx, ty, tz);
                    if (ix<nx-1) tx += .5*hx(ix) + .5*hx(ix+1);
                    ++ip;
                }
                if (iy<ny-1) ty += .5*hy(iy) + .5*hy(iy+1);
            }
            if (iz < nz) tz += hz(iz);
        }

        return ;
    }		// -----  end of method EMSchur3DBase::FillPoints  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  FillSourceTerms
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::FillSourceTerms ( Eigen::Ref<VectorXcr>  ms,
                                          Eigen::Ref<VectorXcr>  Se, Eigen::Ref<VectorXcr> E0,
                                          std::shared_ptr<FieldPoints> dpoint, const Real& omega ) {

        //Complex hsig(0);
        //Complex hsigp(0);
        Real hsig(0);
        Real hsigp(0);
        Real EPS(0); //1e-24);
        Real EPS2(1e-18);
        int iv = 0;
        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny;   ++iy) {
        for (int ix=0; ix<nx+1; ++ix) {

            // Harmonically average sigma
            if (ix == 0) {
                hsig = sigma[ix][iy][iz];
                hsigp = sigmap[ix][iy][iz];
            } else if (ix == nx) {
                hsig = sigma[ix-1][iy][iz];
                hsigp = sigmap[ix-1][iy][iz];
            } else {
                hsig = ((hx[ix]+hx[ix-1])/2.) *
                        (1. / ( (hx[ix  ] / (2.*sigma[ix  ][iy][iz] + EPS))  +
                              (  hx[ix-1] / (2.*sigma[ix-1][iy][iz] + EPS))  ) );
                hsigp = ((hx[ix]+hx[ix-1])/2.) *
                        (1. / ( (hx[ix  ] / (2.*sigmap[ix  ][iy][iz] + EPS))  +
                              (  hx[ix-1] / (2.*sigmap[ix-1][iy][iz] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            if (std::abs(hsigp) < EPS2) hsigp = 0;
            //hsig += Complex(0, omega*EPSILON0);
            //ioms(iv) = Complex(0, omega*MU0*hsig);
            //Se(iv) = ioms(iv)*dpoint->GetEfield(0,iv)(0);
            ms(iv) = MU0*hsig;
            Se(iv) = MU0*hsigp*
                     dpoint->GetEfield(0,iv)(0);
            E0(iv) = dpoint->GetEfield(0,iv)(0);
            ++iv;
        }
        }
        }

        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny+1; ++iy) {
        for (int ix=0; ix<nx;   ++ix) {

            // Harmonically average sigma
            if (iy == 0) {
                hsig = sigma[ix][iy][iz];
                hsigp = sigmap[ix][iy][iz];
            } else if (iy == ny) {
                hsig = sigma[ix][iy-1][iz];
                hsigp = sigmap[ix][iy-1][iz];
            } else {
                hsig = ((hy[iy]+hy[iy-1])/2.) *
                        (1. / ( (hy[iy  ] / (2.*sigma[ix][iy  ][iz] + EPS))  +
                              (  hy[iy-1] / (2.*sigma[ix][iy-1][iz] + EPS))  ) );
                hsigp = ((hy[iy]+hy[iy-1])/2.) *
                        (1. / ( (hy[iy  ] / (2.*sigmap[ix][iy  ][iz] + EPS))  +
                              (  hy[iy-1] / (2.*sigmap[ix][iy-1][iz] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            if (std::abs(hsigp) < EPS2) hsigp = 0;
            //hsig += Complex(0, omega*EPSILON0);
            //ioms(iv) = Complex(0, omega*MU0*hsig);
            //Se(iv) = ioms(iv) * dpoint->GetEfield(0,iv)(1);
            ms(iv) = MU0*hsig;
            Se(iv) = MU0*hsigp* // TODO was hsigp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                     dpoint->GetEfield(0,iv)(1);
            E0(iv) = dpoint->GetEfield(0,iv)(1);
            ++iv;
        }
        }
        }

        for (int iz=0; iz<nz+1; ++iz) {
        for (int iy=0; iy<ny;   ++iy) {
        for (int ix=0; ix<nx;   ++ix) {

            // Harmonically average sigma
            if (iz==0) {
                hsig = sigma[ix][iy][iz];
                hsigp = sigma[ix][iy][iz];
            } else if (iz == nz) {
                hsig = sigma[ix][iy][nz-1];
                hsigp = sigmap[ix][iy][nz-1];
            } else {
                hsig = ((hz[iz]+hz[iz-1])/2.) *
                        (1. / ( (hz[iz  ] / (2.*sigma[ix][iy][iz  ] + EPS))  +
                              (  hz[iz-1] / (2.*sigma[ix][iy][iz-1] + EPS))  ) );
                hsigp = ((hz[iz]+hz[iz-1])/2.) *
                        (1. / ( (hz[iz  ] / (2.*sigmap[ix][iy][iz  ] + EPS))  +
                              (  hz[iz-1] / (2.*sigmap[ix][iy][iz-1] + EPS))  ) );
            }
            if (std::abs(hsig) < EPS2) hsig = 0;
            if (std::abs(hsigp) < EPS2) hsigp = 0;
            //hsig += Complex(0, omega*EPSILON0);
            //ioms(iv) = Complex(0, omega*MU0*hsig);
            //Se(iv) = ioms(iv) * dpoint->GetEfield(0, iv)(2);
            ms(iv) = MU0*hsig;
            Se(iv) = MU0*hsig
                     *dpoint->GetEfield(0, iv)(2);
            E0(iv) = dpoint->GetEfield(0, iv)(2);
            ++iv;
        }
        }
        }

        return ;
    }		// -----  end of method EMSchur3DBase::FillSourceTerms  -----

    VectorXcr EMSchur3DBase::StaggeredGridCurl( Eigen::Ref<VectorXcr> A ) {
        VectorXcr B    = VectorXcr::Zero(3*nx*ny*nz);
        VectorXcr dzdy = VectorXcr::Zero(  nx*ny*nz);
        VectorXcr dydz = VectorXcr::Zero(  nx*ny*nz);
        VectorXcr dxdz = VectorXcr::Zero(  nx*ny*nz);
        VectorXcr dzdx = VectorXcr::Zero(  nx*ny*nz);
        VectorXcr dydx = VectorXcr::Zero(  nx*ny*nz);
        VectorXcr dxdy = VectorXcr::Zero(  nx*ny*nz);
        /*
              Curl_x = dzdy - dydz
              Curl_y = dxdz - dzdx
              Curl_z = dydx - dxdy
         */
        ////////////////////////////////
        // X component
        int ii = 0;
        int iix = 0;
        for (int iz=0; iz<nz; ++iz) {
        for (int iy=0; iy<ny; ++iy) {
            for (int ix=0; ix<nx; ++ix) {
                dxdz(iix)  = ( A(ii+(nx+1)*ny) - A(ii) ) / hz[iz];
                dxdy(iix)  = ( A(ii+(nx+1)   ) - A(ii) ) / hy[iy];
                ++ii;
                ++iix;
            }
            // Jump around boundary
            ++ii;
        }
        }

        ///////////////////////////////
        // Y component
        int iiy = 0;
        for (int iz=0; iz<nz;   ++iz) {
            for (int iy=0; iy<ny; ++iy) {
                for (int ix=0; ix<nx;   ++ix) {
                    dydz(iiy)  = ( A(ii+(nx)*(ny+1)) - A(ii) ) / hz[iz];
                    dydx(iiy)  = ( A(ii+1          ) - A(ii) ) / hx[ix];
                    ++iiy;
                    ++ii;
                }
            }
            // Jump around boundary
            ii += nx;
        }

        ////////////////////////////
        // Z component
        int iiz = 0;
        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny; ++iy) {
            for (int ix=0; ix<nx;   ++ix) {
                //Az(iiz)    = 0.5 * ( A(ii)    + A(ii+nx*ny));
                dzdy(iiz)  = ( A(ii+nx       ) - A(ii) ) / hy[iy];
                dzdx(iiz)  = ( A(ii+1        ) - A(ii) ) / hx[ix];
                ++iiz;
                ++ii;
            }
            //++iiz;
            //ii += nx*ny;
        }
        }

        B.segment(         0, nx*ny*nz) = dzdy - dydz;
        B.segment(  nx*ny*nz, nx*ny*nz) = dxdz - dzdx;
        B.segment(2*nx*ny*nz, nx*ny*nz) = dydx - dxdy;

        return B;
    }

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  WriteVTKResults
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::WriteVTKResults ( const std::string& fname, Eigen::Ref<VectorXcr> A,
            Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0, Eigen::Ref<VectorXcr> E,
            Eigen::Ref<VectorXcr> Phi, Eigen::Ref<VectorXcr> ADiv, Eigen::Ref<VectorXcr> ADiv2,
            Eigen::Ref<VectorXcr> B ) {

        auto VTKGridExporter = RectilinearGridVTKExporter::NewSP();
        VTKGridExporter->SetGrid(Grid);

        // VTK arrays to fill
        // scalars
        vtkDoubleArray *sigmaArray = vtkDoubleArray::New();
        vtkDoubleArray *sigmapArray = vtkDoubleArray::New();
        vtkDoubleArray *phiRealArray = vtkDoubleArray::New();
        vtkDoubleArray *phiImagArray = vtkDoubleArray::New();
        vtkDoubleArray *ADivRealArray = vtkDoubleArray::New();
        vtkDoubleArray *ADivImagArray = vtkDoubleArray::New();
        vtkDoubleArray *ADiv2RealArray = vtkDoubleArray::New();
        vtkDoubleArray *ADiv2ImagArray = vtkDoubleArray::New();

        // vectors
        vtkDoubleArray *AReal = vtkDoubleArray::New();
        vtkDoubleArray *AImag = vtkDoubleArray::New();
        vtkDoubleArray *BReal = vtkDoubleArray::New();
        vtkDoubleArray *BImag = vtkDoubleArray::New();
        vtkDoubleArray *SeReal = vtkDoubleArray::New();
        vtkDoubleArray *SeImag = vtkDoubleArray::New();
        vtkDoubleArray *E0Real = vtkDoubleArray::New();
        vtkDoubleArray *E0Imag = vtkDoubleArray::New();
        vtkDoubleArray *EReal = vtkDoubleArray::New();
        vtkDoubleArray *EImag = vtkDoubleArray::New();

        AReal->SetNumberOfComponents(3);
        AImag->SetNumberOfComponents(3);
        BReal->SetNumberOfComponents(3);
        BImag->SetNumberOfComponents(3);
        SeReal->SetNumberOfComponents(3);
        SeImag->SetNumberOfComponents(3);
        E0Real->SetNumberOfComponents(3);
        E0Imag->SetNumberOfComponents(3);
        EReal->SetNumberOfComponents(3);
        EImag->SetNumberOfComponents(3);

        // Interpolate
        Eigen::VectorXcd  Ax(nx*ny*nz);       // A
        Eigen::VectorXcd  Ay(nx*ny*nz);
        Eigen::VectorXcd  Az(nx*ny*nz);
        Eigen::VectorXcd  Bx(nx*ny*nz);       // B
        Eigen::VectorXcd  By(nx*ny*nz);
        Eigen::VectorXcd  Bz(nx*ny*nz);
        Eigen::VectorXcd  Sex(nx*ny*nz);      // Se
        Eigen::VectorXcd  Sey(nx*ny*nz);
        Eigen::VectorXcd  Sez(nx*ny*nz);
        Eigen::VectorXcd  E0x(nx*ny*nz);      // E0
        Eigen::VectorXcd  E0y(nx*ny*nz);
        Eigen::VectorXcd  E0z(nx*ny*nz);
        Eigen::VectorXcd  Ex(nx*ny*nz);       // E
        Eigen::VectorXcd  Ey(nx*ny*nz);
        Eigen::VectorXcd  Ez(nx*ny*nz);

        ////////////////////////////////
        // X component
        int ii = 0;
        int iix = 0;
        int ic = 0;
        for (int iz=0; iz<nz; ++iz) {
        for (int iy=0; iy<ny; ++iy) {
            for (int ix=0; ix<nx; ++ix) {
                Ax(iix)    = 0.5 * ( A(ii) +  A(ii+1));
                Sex(iix)   = 0.5 * (Se(ii) + Se(ii+1));
                E0x(iix)   = 0.5 * (E0(ii) + E0(ii+1));
                Ex(iix)    = 0.5 * ( E(ii) +  E(ii+1));
                Bx(iix)    = B(ic);
                ++ii;
                ++iix;
                ++ic;
            }
            // Jump around boundary
            ++ii;
        }
        }

        ///////////////////////////////
        // Y component
        int iiy = 0;
        for (int iz=0; iz<nz;   ++iz) {
            for (int iy=0; iy<ny; ++iy) {
                for (int ix=0; ix<nx;   ++ix) {
                    Ay(iiy)    = 0.5 * ( A(ii) +  A(ii+nx));
                    Sey(iiy)   = 0.5 * (Se(ii) + Se(ii+nx));
                    E0y(iiy)   = 0.5 * (E0(ii) + E0(ii+nx));
                    Ey(iiy)    = 0.5 * ( E(ii) +  E(ii+nx));
                    By(iiy)    = B(ic);
                    ++iiy;
                    ++ii;
                    ++ic;
                }
            }
            // Jump around boundary
            ii += nx;
        }

        ////////////////////////////
        // Z component
        int iiz = 0;
        for (int iz=0; iz<nz;   ++iz) {
        for (int iy=0; iy<ny; ++iy) {
            for (int ix=0; ix<nx;   ++ix) {
                Az(iiz)    = 0.5 * ( A(ii)    + A(ii+nx*ny));
                Sez(iiz)   = 0.5 * (Se(ii)    + Se(ii+nx*ny));
                E0z(iiz)   = 0.5 * (E0(ii)    + E0(ii+nx*ny));
                Ez(iiz)    = 0.5 * ( E(ii)    + E(ii+nx*ny));
                Bz(iiz)    = B(ic);
                ++iiz;
                ++ii;
                ++ic;
            }
            //++iiz;
            //ii += nx*ny;
        }
        }

        //VTKGridExporter->GetVTKGrid();
        int i = 0;
        for (int iz=0; iz<nz; ++iz) {
        for (int iy=0; iy<ny; ++iy) {
        for (int ix=0; ix<nx; ++ix) {
            // scalars
            //sigmaArray->InsertTuple1(i,  std::real(sigma[ix][iy][iz]) );
            //sigmapArray->InsertTuple1(i,  std::real(sigmap[ix][iy][iz]) );
            sigmaArray->InsertTuple1(i,  sigma[ix][iy][iz] );
            sigmapArray->InsertTuple1(i, sigmap[ix][iy][iz] );

            phiRealArray->InsertTuple1(i, std::real(Phi(i)));
            phiImagArray->InsertTuple1(i, std::imag(Phi(i)));
            ADivRealArray->InsertTuple1(i, std::real(ADiv(i)));
            ADivImagArray->InsertTuple1(i, std::imag(ADiv(i)));
            ADiv2RealArray->InsertTuple1(i, std::real(ADiv2(i)));
            ADiv2ImagArray->InsertTuple1(i, std::imag(ADiv2(i)));

            // vectors
            AReal->InsertTuple3(i, std::real(Ax(i)), std::real(Ay(i)), std::real(Az(i)));
            AImag->InsertTuple3(i, std::imag(Ax(i)), std::imag(Ay(i)), std::imag(Az(i)));
            BReal->InsertTuple3(i, std::real(Bx(i)), std::real(By(i)), std::real(Bz(i)));
            BImag->InsertTuple3(i, std::imag(Bx(i)), std::imag(By(i)), std::imag(Bz(i)));
            SeReal->InsertTuple3(i, std::real(Sex(i)), std::real(Sey(i)), std::real(Sez(i)));
            SeImag->InsertTuple3(i, std::imag(Sex(i)), std::imag(Sey(i)), std::imag(Sez(i)));
            E0Real->InsertTuple3(i, std::real(E0x(i)), std::real(E0y(i)), std::real(E0z(i)));
            E0Imag->InsertTuple3(i, std::imag(E0x(i)), std::imag(E0y(i)), std::imag(E0z(i)));
            EReal-> InsertTuple3(i, std::real(Ex(i)), std::real(Ey(i)), std::real(Ez(i)));
            EImag-> InsertTuple3(i, std::imag(Ex(i)), std::imag(Ey(i)), std::imag(Ez(i)));

            ++i;
        }
        }
        }

        AReal->SetName("A_real");
        AImag->SetName("A_imag");
        BReal->SetName("B_real");
        BImag->SetName("B_imag");
        SeReal->SetName("Se_real");
        SeImag->SetName("Se_imag");
        E0Real->SetName("E0_real");
        E0Imag->SetName("E0_imag");
        EReal->SetName("E_real");
        EImag->SetName("E_imag");

        phiRealArray->SetName("phi Real");
        phiImagArray->SetName("phi imag");
        ADivRealArray->SetName("ini div A real");
        ADivImagArray->SetName("ini div A imag");
        ADiv2RealArray->SetName("div A real");
        ADiv2ImagArray->SetName("div A imag");

        // Add scalar fields
        sigmaArray->SetName("sigma");
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(sigmaArray);

        sigmapArray->SetName("sigma prime");
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(sigmapArray);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(phiRealArray);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(phiImagArray);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADivRealArray);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADivImagArray);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADiv2RealArray);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(ADiv2ImagArray);

        // Add Vector Arrays
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(AReal);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(AImag);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(SeReal);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(SeImag);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(E0Real);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(E0Imag);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(EReal);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(EImag);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(BReal);
        VTKGridExporter->GetVTKGrid()->GetCellData()->AddArray(BImag);

        // Write

        vtkXMLRectilinearGridWriter *gridWrite = vtkXMLRectilinearGridWriter::New();
            gridWrite->SetInputData( VTKGridExporter->GetVTKGrid() );
            gridWrite->SetFileName( (fname + std::string(".vtr")).c_str());
            gridWrite->Update();
            gridWrite->Write();

        sigmaArray->Delete();
        gridWrite->Delete();

        phiRealArray->Delete();
        phiImagArray->Delete();
        ADivRealArray->Delete();
        ADivImagArray->Delete();
        ADiv2RealArray->Delete();
        ADiv2ImagArray->Delete();

        // TODO delete vectors too!
        AReal->Delete();
        AImag->Delete();
        BReal->Delete();
        BImag->Delete();
        SeReal->Delete();
        SeImag->Delete();
        E0Real->Delete();
        E0Imag->Delete();
        EReal->Delete();
        EImag->Delete();

        return ;
    }		// -----  end of method EMSchur3DBase::WriteVTKResults  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  LoadMeshToolsConductivityModel
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::LoadMeshToolsConductivityModel ( const std::string& fname ) {

        if (sigma) Delete3DScalar(sigma);
        //Allocate3DScalar(sigma, Complex(0));
        Allocate3DScalar(sigma, 0.);

        // TODO, there are smarter ways to do this, since the bacground is 1D. But for now, we just have
        //       them the same. Eases logic in harmonic averaging. And it's just in one thread, so not so
        //       hateful
        if (!LayModel) {
            std::cerr << "In EMSchur3DBase::LoadMeshToolsConductivityModel, you need to set your bacground 1D model\n";
            exit(EXIT_FAILURE);
        }

        if (sigmap) Delete3DScalar(sigmap);
        //Allocate3DScalar(sigmap, Complex(0));
        Allocate3DScalar(sigmap, 0.);


        auto Parser = ASCIIParser::NewSP();
        Parser->Open(fname);
        std::vector<Real> sigmod = Parser->ReadReals( -1 );
        if ( static_cast<int>(sigmod.size()) != nx*ny*nz) {
            std::cerr << "In EMSchur3DBase::LoadMeshToolsConductivityModel model sizes are not in agreement\n";
            exit(EXIT_FAILURE);
        }
        Parser->Close();


        int ic(0);
        for (int ix=0; ix<nx; ++ix) {
        for (int iy=0; iy<ny; ++iy) {
        Real pz = Grid->GetOz();
        for (int iz=0; iz<nz; ++iz) {
            if ( sigmod[ic] != sigmod[ic] ) {
                std::cout << "sigmod problems " << ic << std::endl;
            }
            sigma[ix][iy][iz]  = sigmod[ic];
            //sigmap[ix][iy][iz] = sigmod[ic] - LayModel->GetLayerConductivity(LayModel->GetLayerAtThisDepth(pz));
            sigmap[ix][iy][iz] = sigmod[ic] - LayModel->GetLayerConductivity(LayModel->GetLayerAtThisDepth(pz)).real() ;

            ++ic;
            pz += hz[iz];
        }
        }
        }

        // Marker and cell
        MAC  = VectorXi::Zero(nx*ny*nz);
        idx.clear();
        ic = 0;
        for (int iz=0; iz<nz; ++iz) {
        for (int iy=0; iy<ny; ++iy) {
        for (int ix=0; ix<nx; ++ix) {
            if (sigma[ix][iy][iz] > 0.) {
                MAC(ic) = 1;
                idx.push_back(ic);
            }
            ++ic;
        }
        }
        }

        return ;
    }

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  SetAEMSurvey
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::SetAEMSurvey ( std::shared_ptr<AEMSurvey> SurveyIn ) {
        Survey = SurveyIn;
        // determine how many frequencies we are dealing with
        Omegas = 2.*PI*Survey->GetFrequencies();
        return ;
    }		// -----  end of method EMSchur3DBase::SetAEMSurvey  -----// -----  end of method EMSchur3DBase::LoadMeshToolsConductivityModel  -----

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3DBase
    //      Method:  Solve
    //--------------------------------------------------------------------------------------
    void EMSchur3DBase::Solve (  ) {

        // TODO, these should throw exceptions for a GUI to catch
        if (!LayModel or !Survey or !Grid or !sigma) {
            std::cerr << "You need to set something. EMSChur3D::Solve()";
            exit(EXIT_FAILURE);
        }

        //BuildD(); // strangely necessary, bug, track down.
        Setup();

        std::cout << "Entering parallel solve" << std::endl;
        //#ifdef LEMMAUSEOMP
        //#pragma omp parallel for schedule(static,1)
        //#endif
        for (int isource=0; isource<Survey->GetNumberOfSources(); ++isource) {
            SolveSource(Survey->GetSource(isource), isource);
        }

        return ;
    }		// -----  end of method EMSchur3DBase::Solve  -----

}		// -----  end of Lemma  name  -----