Added coupling routine

examples/CMakeLists.txt

@@ -1,11 +1,16 @@
 add_executable( KernelV0 KernelV0.cpp  )
 target_link_libraries(  KernelV0  "lemmacore" "fdem1d" "merlin")
 
+add_executable( Coupling Coupling.cpp  )
+target_link_libraries(  Coupling  "lemmacore" "fdem1d" "merlin")
+
 # Linking
 if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT ) 
 	target_link_libraries( KernelV0 ${VTK_LIBRARIES})
+	target_link_libraries( Coupling ${VTK_LIBRARIES})
 endif()
 
 INSTALL_TARGETS( "/share/Merlin/"
 	KernelV0
+	Coupling
 )

examples/Coupling.cpp

@@ -0,0 +1,101 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      11/11/2016 02:44:37 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Lemma Software, LLC
+ */
+
+#include <Merlin>
+using namespace Lemma;
+
+std::shared_ptr<PolygonalWireAntenna> CircularLoop ( int nd, Real radius, Real Offsetx, Real Offsety ) ;
+void MoveLoop( std::shared_ptr<PolygonalWireAntenna> Loop, int nd, Real Radius, Real Offsetx, Real Offsety );
+
+int main(int argc, char** argv) {
+
+    Real offset = atof(argv[1]);
+
+	auto earth = LayeredEarthEM::NewSP();
+		earth->SetNumberOfLayers(3);
+		earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
+		earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
+        // Set mag field info
+        // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
+        earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
+
+    // Transmitter loops
+    auto Tx1 = CircularLoop(21, 15, 100, 100);
+    auto Tx2 = CircularLoop(21, 15, 100, 100 + offset);
+
+    auto Kern = Coupling::NewSP();
+        Kern->PushCoil( "Coil 1", Tx1 );
+        Kern->PushCoil( "Coil 2", Tx2 );
+        Kern->SetLayeredEarthEM( earth );
+
+        Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
+        Kern->SetIntegrationOrigin( (Vector3r() << 0,0,2).finished() );
+        Kern->SetTolerance( 1e-7 ); // 1e-12
+
+    std::vector<std::string> tx = {std::string("Coil 1")};
+    std::vector<std::string> rx = {std::string("Coil 2")};
+    VectorXr Offsets = VectorXr::LinSpaced(60,0.25,60); // nbins, low, high
+
+    auto outfile = std::ofstream("coupling.dat");
+    for (int ii=0; ii< Offsets.size(); ++ii) {
+        MoveLoop(Tx2, 21, 15, 100, 100 + Offsets(ii));
+        Complex coupling = Kern->Calculate( tx, rx, false );
+        std::cout << "coupling " << coupling << std::endl;
+        outfile << Offsets(ii) << "\t" <<  std::real(coupling) << "\t" << std::imag(coupling) << std::endl;
+    }
+    outfile.close();
+}
+
+std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {
+
+    auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
+         Tx1->SetNumberOfPoints(nd);
+
+    VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
+    int ii;
+    for (ii=0; ii<nd; ++ii) {
+        Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)),  -1e-3));
+    }
+    //Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety,  -1e-3));
+
+    Tx1->SetCurrent(1.);
+    Tx1->SetNumberOfTurns(1);
+    Tx1->SetNumberOfFrequencies(1);
+    Tx1->SetFrequency(0,2246);
+
+    return Tx1;
+}
+
+void MoveLoop( std::shared_ptr<Lemma::PolygonalWireAntenna> Tx1, int nd, Real Radius, Real Offsetx, Real Offsety ) {
+
+    Tx1->SetNumberOfPoints(nd);
+
+    VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
+    int ii;
+    for (ii=0; ii<nd; ++ii) {
+        Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)),  -1e-3));
+    }
+    //Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety,  -1e-3));
+
+    Tx1->SetCurrent(1.);
+    Tx1->SetNumberOfTurns(1);
+    Tx1->SetNumberOfFrequencies(1);
+    Tx1->SetFrequency(0,2246);
+
+}

examples/KernelV0.cpp

@@ -47,10 +47,12 @@ int main(int argc, char** argv) {
 
         Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
         Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
-        Kern->SetTolerance( 1e-12 );
+        Kern->SetTolerance( 1e-12 ); // 1e-12
 
         Kern->SetPulseDuration(0.020);
         VectorXr I(36);
+
+        // off from VC by 1.075926340216996
         // Pulses from Wyoming Red Buttes exp 0
         I << 397.4208916184016, 352.364477036168, 313.0112765842783, 278.37896394065376, 247.81424224324982,
              220.77925043190442, 196.76493264105017, 175.31662279234038, 156.0044839325404, 138.73983004230124,
@@ -63,22 +65,22 @@ int main(int argc, char** argv) {
         Kern->SetPulseCurrent( I ); // nbins, low, high
 
         //Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 30, 3, 45.5 ) ); // nlay, low, high
-        VectorXr interfaces = VectorXr::LinSpaced( 31, 1, 45.5 ); // nlay, low, high
+        VectorXr interfaces = VectorXr::LinSpaced( 41, .5, 45.5 ); // nlay, low, high
         Real thick = .5;
         for (int ilay=1; ilay<interfaces.size(); ++ilay) {
             interfaces(ilay) = interfaces(ilay-1) + thick;
-            thick *= 1.1;
+            thick *= 1.05;
         }
         Kern->SetDepthLayerInterfaces( interfaces ); // nlay, low, high
 
     // We could, I suppose, take the earth model in here? For non-linear that
     // may be more natural to work with?
     std::vector<std::string> tx = {std::string("Coil 1"), std::string("Coil 2") };
-    //std::vector<std::string> tx = {std::string("Coil 1")};
     std::vector<std::string> rx = {std::string("Coil 1")};
     Kern->CalculateK0( tx, rx, false );
 
     std::ofstream dout = std::ofstream(std::string("k-")+ std::string(argv[1])+ std::string(".dat"));
+    //std::ofstream dout = std::ofstream(std::string("k-coincident.dat"));
         dout << interfaces.transpose() << std::endl;
         dout << I.transpose() << std::endl;
         dout << "#real\n";
@@ -88,6 +90,7 @@ int main(int argc, char** argv) {
         dout.close();
 
     std::ofstream out = std::ofstream(std::string("k-")+std::string(argv[1])+std::string(".yaml"));
+    //std::ofstream out = std::ofstream(std::string("k-coincident.yaml"));
     out << *Kern;
     out.close();
 }

examples/plotkernel.py

@@ -22,7 +22,7 @@ ax1 =  fig.add_axes( [.100,.125,.335,.775] )
 ax2 =  fig.add_axes( [.465,.125,.335,.775] , sharex=ax1, sharey=ax1 )
 axcb = fig.add_axes( [.835,.125,.040,.775] )
 
-ccmap = "seismic" # RdBu
+ccmap = "coolwarm" # RdBu
 
 # Real plot
 ax1.pcolormesh(X, Y, K[0:nx//2].T, cmap=ccmap, vmin=-np.max(np.abs(K)), vmax=np.max(np.abs(K)))

include/Coupling.h

@@ -0,0 +1,258 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      11/11/2016 01:47:34 PM
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Lemma Software, LLC
+ * @copyright Copyright (c) 2008, Colorado School of Mines
+ */
+#ifndef  COUPLING_INC
+#define  COUPLING_INC
+
+#pragma once
+#include "LemmaObject.h"
+#include "LayeredEarthEM.h"
+#include "PolygonalWireAntenna.h"
+#include "EMEarth1D.h"
+
+#ifdef LEMMAUSEVTK
+#include "vtkHyperOctree.h"
+#include "vtkHyperOctreeCursor.h"
+#include "vtkXMLHyperOctreeWriter.h"
+#include "vtkDoubleArray.h"
+#endif
+
+namespace Lemma {
+
+    /**
+     * \ingroup Merlin
+     * \brief
+     * \details
+     */
+    class Coupling : public LemmaObject {
+
+        friend std::ostream &operator<<(std::ostream &stream, const Coupling &ob);
+
+        protected:
+        /*
+         *  This key is used to lock the constructor. It is protected so that inhereted
+         *  classes also have the key to contruct their base class.
+         */
+        struct ctor_key {};
+
+        public:
+
+        // ====================  LIFECYCLE     =======================
+
+        /**
+         * Default constructor.
+         * @note This method is locked, and cannot be called directly.
+         *       The reason that the method is public is to enable the use
+         *       of make_shared whilst enforcing the use of shared_ptr,
+         *       in c++-17, this curiosity may be resolved.
+         * @see Coupling::NewSP
+         */
+        explicit Coupling ( const ctor_key& );
+
+        /**
+         * DeSerializing constructor.
+         * @note This method is locked, and cannot be called directly.
+         *       The reason that the method is public is to enable the use
+         *       of make_shared whilst enforcing the use of shared_ptr,
+         *       in c++-17, this curiosity may be resolved.
+         * @see Coupling::DeSerialize
+         */
+        Coupling ( const YAML::Node& node, const ctor_key& );
+
+        /**
+         * Default destructor.
+         * @note This method should never be called due to the mandated
+         *       use of smart pointers. It is necessary to keep the method
+         *       public in order to allow for the use of the more efficient
+         *       make_shared constructor.
+         */
+        virtual ~Coupling ();
+
+        /**
+         *  Uses YAML to serialize this object.
+         *  @return a YAML::Node
+         *  @see Coupling::DeSerialize
+         */
+        virtual YAML::Node Serialize() const;
+
+        /*
+         *  Factory method for generating concrete class.
+         *  @return a std::shared_ptr of type Coupling
+         */
+        static std::shared_ptr< Coupling > NewSP();
+
+        /**
+         *   Constructs an Coupling object from a YAML::Node.
+         *   @see Coupling::Serialize
+         */
+        static std::shared_ptr<Coupling> DeSerialize(const YAML::Node& node);
+
+        // ====================  OPERATORS     =======================
+
+        // ====================  OPERATIONS    =======================
+
+        /**
+         * @return std::shared_ptr<LayeredEarthEM>
+         */
+        inline std::shared_ptr<LayeredEarthEM> GetSigmaModel (  ) {
+            return SigmaModel;
+        }		// -----  end of method Coupling::get_SigmaModel  -----
+
+        /**
+         * @param[in] value the 1D-EM model used for calculations
+         */
+        inline void SetLayeredEarthEM ( std::shared_ptr< LayeredEarthEM > value ) {
+            SigmaModel	= value;
+            return ;
+        }		// -----  end of method Coupling::set_SigmaModel  -----
+
+        /**
+         *
+         */
+        inline void SetIntegrationSize ( const Vector3r& size ) {
+            Size = size;
+            return ;
+        }		// -----  end of method Coupling::SetIntegrationSize  -----
+
+        /**
+         *
+         */
+        inline void SetIntegrationOrigin ( const Vector3r& origin ) {
+            Origin = origin;
+            return ;
+        }		// -----  end of method Coupling::SetIntegrationOrigin  -----
+
+        /**
+         *   Assign transmiter coils
+         */
+        inline void PushCoil( const std::string& label, std::shared_ptr<PolygonalWireAntenna> ant ) {
+            TxRx[label] = ant;
+        }
+
+        /**
+         *   Calculates a single imaging kernel, however, phased arrays are supported
+         *   so that more than one transmitter and/or receiver can be specified.
+         *   @param[in] tx is the list of transmitters to use for a kernel, use the same labels as
+         *              used in PushCoil.
+         *   @param[in] rx is the list of receivers to use for a kernel, use the same labels as
+         *              used in PushCoil. @see PushCoil
+         *   @param[in] vtkOutput generates a VTK hyperoctree file as well, useful for visualization.
+         *              requires compilation of Lemma with VTK.
+         */
+        Complex Calculate (const std::vector< std::string >& tx, const std::vector< std::string >& rx,
+                bool vtkOutput=false );
+
+        /**
+         *  Sets the tolerance to use for making the adaptive mesh
+         *
+         */
+        inline void SetTolerance(const Real& ttol) {
+            tol = ttol;
+        }
+
+        inline void SetPulseDuration(const Real& taup) {
+            Taup = taup;
+        }
+
+        // ====================  INQUIRY       =======================
+        /**
+         *  Returns the name of the underlying class, similiar to Python's type
+         *  @return string of class name
+         */
+        virtual inline std::string GetName() const {
+            return CName;
+        }
+
+        protected:
+
+        // ====================  LIFECYCLE     =======================
+
+        /** Copy is disabled */
+        Coupling( const Coupling& ) = delete;
+
+        private:
+
+        /**
+         *  Returns the kernel value for an input prism
+         */
+        Complex f( const Vector3r& r, const Real& volume , const Vector3cr& Ht, const Vector3cr& Hr);
+
+        void IntegrateOnOctreeGrid( bool vtkOutput=false );
+
+        /**
+         *  Recursive call to integrate a function on an adaptive Octree Grid.
+         *  For efficiency's sake the octree grid is not stored, as only the
+         *  integral (sum) is of interest. The logic for grid refinement is based
+         *  on an Octree representation of the domain. If an Octree representation
+         *  of the kernel is desired, call alternative version @see EvaluateKids2
+         *  @param[in] size gives the domain size, in  metres
+         *  @param[in] level gives the current level of the octree grid, call with 0 initially
+         *  @param[in] cpos is the centre position of the parent cuboid
+         */
+        void EvaluateKids(  const Vector3r& size, const int& level, const Vector3r& cpos,
+                            const Complex& parentVal );
+
+        #ifdef LEMMAUSEVTK
+        /**
+         *  Same functionality as @see EvaluateKids, but includes generation of a VTK
+         *  HyperOctree, which is useful for visualization.
+         */
+        void EvaluateKids2(  const Vector3r& size, const int& level, const Vector3r& cpos,
+                            const VectorXcr& parentVal, vtkHyperOctree* octree, vtkHyperOctreeCursor* curse );
+
+        void GetPosition( vtkHyperOctreeCursor* Cursor, Real* p );
+        #endif
+
+        // ====================  DATA MEMBERS  =========================
+
+        int                                       ilay;
+        int                                       nleaves;
+        int                                       minLevel=4;
+        int                                       maxLevel=8;
+
+        Real                                      VOLSUM;
+        Real                                      tol=1e-11;
+        Real                                      Taup = .020;  // Sec
+
+        Complex                                   SUM;
+
+        Vector3r                                  Size;
+        Vector3r                                  Origin;
+
+        std::shared_ptr< LayeredEarthEM >         SigmaModel = nullptr;
+        std::shared_ptr< FieldPoints >            cpoints;
+
+        std::map< std::string , std::shared_ptr< PolygonalWireAntenna > >  TxRx;
+        std::map< std::string , std::shared_ptr< EMEarth1D > >             EMEarths;
+
+        #ifdef LEMMAUSEVTK
+        std::map< int, VectorXcr  >               LeafDict;
+        std::map< int, int     >                  LeafDictIdx;
+        std::map< int, Real     >                 LeafDictErr;
+        #endif
+
+        /** ASCII string representation of the class name */
+        static constexpr auto CName = "Coupling";
+
+    }; // -----  end of class  Coupling  -----
+}  // -----  end of namespace Lemma ----
+
+/* vim: set tabstop=4 expandtab */
+/* vim: set filetype=cpp */
+
+#endif   // ----- #ifndef COUPLING_INC  -----

include/Merlin

@@ -1,4 +1,5 @@
 #include "KernelV0.h"
+#include "Coupling.h"
 
 /* vim: set tabstop=4 expandtab: */
 /* vim: set filetype=cpp: */

src/CMakeLists.txt

@@ -1,4 +1,5 @@
 set (MERLINSOURCE
 	${CMAKE_CURRENT_SOURCE_DIR}/KernelV0.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/Coupling.cpp
 	PARENT_SCOPE
 )

src/Coupling.cpp

@@ -0,0 +1,444 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      11/11/2016 01:47:25 PM
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Lemma Software, LLC
+ * @copyright Copyright (c) 2008, Colorado School of Mines
+ */
+
+
+#include "Coupling.h"
+#include "FieldPoints.h"
+
+namespace Lemma {
+
+    // ====================  FRIEND METHODS  =====================
+
+    std::ostream &operator << (std::ostream &stream, const Coupling &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc ---
+        return stream;
+    }
+
+    // ====================  LIFECYCLE     =======================
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  Coupling
+    // Description:  constructor (locked)
+    //--------------------------------------------------------------------------------------
+    Coupling::Coupling (const ctor_key&) : LemmaObject( ) {
+
+    }  // -----  end of method Coupling::Coupling  (constructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  Coupling
+    // Description:  DeSerializing constructor (locked)
+    //--------------------------------------------------------------------------------------
+    Coupling::Coupling (const YAML::Node& node, const ctor_key&) : LemmaObject(node) {
+
+    }  // -----  end of method Coupling::Coupling  (constructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  NewSP()
+    // Description:  public constructor returing a shared_ptr
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr< Coupling >  Coupling::NewSP() {
+        return std::make_shared< Coupling >( ctor_key() );
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  ~Coupling
+    // Description:  destructor (protected)
+    //--------------------------------------------------------------------------------------
+    Coupling::~Coupling () {
+
+    }  // -----  end of method Coupling::~Coupling  (destructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  Serialize
+    //--------------------------------------------------------------------------------------
+    YAML::Node  Coupling::Serialize (  ) const {
+        YAML::Node node = LemmaObject::Serialize();
+        node.SetTag( GetName() );
+
+        // Coils Transmitters & Receivers
+        for ( auto txm : TxRx) {
+            node[txm.first] = txm.second->Serialize();
+        }
+        // LayeredEarthEM
+        node["SigmaModel"] = SigmaModel->Serialize();
+
+        node["tol"] = tol;
+        node["minLevel"] = minLevel;
+        node["maxLevel"] = maxLevel;
+
+        return node;
+    }		// -----  end of method Coupling::Serialize  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  DeSerialize
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr<Coupling> Coupling::DeSerialize ( const YAML::Node& node  ) {
+        if (node.Tag() !=  "Coupling" ) {
+            throw  DeSerializeTypeMismatch( "Coupling", node.Tag());
+        }
+        return std::make_shared< Coupling > ( node, ctor_key() );
+    }		// -----  end of method Coupling::DeSerialize  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  DeSerialize
+    //--------------------------------------------------------------------------------------
+    Complex Coupling::Calculate (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
+            bool vtkOutput ) {
+        // All EM calculations will share same field points
+        cpoints = FieldPoints::NewSP();
+            cpoints->SetNumberOfPoints(8);
+        for (auto tx : Tx) {
+            // Set up EMEarth
+            EMEarths[tx] = EMEarth1D::NewSP();
+                EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
+                EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
+                EMEarths[tx]->AttachFieldPoints( cpoints );
+         		EMEarths[tx]->SetFieldsToCalculate(H);
+                // TODO query for method, altough with flat antennae, this is fastest
+                EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
+                EMEarths[tx]->SetTxRxMode(TX);
+                TxRx[tx]->SetCurrent(1.);
+        }
+        for (auto rx : Rx) {
+            if (EMEarths.count(rx)) {
+                EMEarths[rx]->SetTxRxMode(TXRX);
+            } else {
+                EMEarths[rx] = EMEarth1D::NewSP();
+                    EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
+                    EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
+                    EMEarths[rx]->AttachFieldPoints( cpoints );
+         		    EMEarths[rx]->SetFieldsToCalculate(H);
+                    // TODO query for method, altough with flat antennae, this is fastest
+                    EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
+                    EMEarths[rx]->SetTxRxMode(RX);
+                    TxRx[rx]->SetCurrent(1.);
+            }
+        }
+        SUM = 0;
+        IntegrateOnOctreeGrid( vtkOutput );
+        std::cout << "\nFinished KERNEL\n";
+        return SUM;
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  IntegrateOnOctreeGrid
+    //--------------------------------------------------------------------------------------
+    void Coupling::IntegrateOnOctreeGrid( bool vtkOutput) {
+
+        Vector3r cpos = Origin + Size/2.;
+
+        VOLSUM = 0;
+        nleaves = 0;
+        if (!vtkOutput) {
+            EvaluateKids( Size, 0, cpos, Complex(100.));
+        } else {
+        #ifdef LEMMAUSEVTK
+            vtkHyperOctree* oct = vtkHyperOctree::New();
+                oct->SetDimension(3);
+                oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
+                oct->SetSize( Size(0), Size(1), Size(2) );
+            vtkHyperOctreeCursor* curse = oct->NewCellCursor();
+                curse->ToRoot();
+            EvaluateKids2( Size, 0, cpos, VectorXcr::Ones(PulseI.size()), oct, curse );
+
+            for (int iq=0; iq<PulseI.size(); ++iq) {
+
+            // Fill in leaf data
+            vtkDoubleArray* kr = vtkDoubleArray::New();
+                kr->SetNumberOfComponents(1);
+                kr->SetName("Re($K_0$)");
+                kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            vtkDoubleArray* ki = vtkDoubleArray::New();
+                ki->SetNumberOfComponents(1);
+                ki->SetName("Im($K_0$)");
+                ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            vtkDoubleArray* km = vtkDoubleArray::New();
+                km->SetNumberOfComponents(1);
+                km->SetName("mod($K_0$)");
+                km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            vtkIntArray* kid = vtkIntArray::New();
+                kid->SetNumberOfComponents(1);
+                kid->SetName("ID");
+                kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            vtkIntArray* kerr = vtkIntArray::New();
+                kerr->SetNumberOfComponents(1);
+                kerr->SetName("nleaf");
+
+            //Real LeafVol(0);
+            for (auto leaf : LeafDict) {
+                kr->InsertTuple1( leaf.first, std::real(leaf.second(iq)) );
+                ki->InsertTuple1( leaf.first, std::imag(leaf.second(iq)) );
+                km->InsertTuple1( leaf.first, std::abs(leaf.second(iq)) );
+                kid->InsertTuple1( leaf.first, leaf.first );
+                //LeafVol += std::real(leaf.second);
+            }
+            //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
+
+            for (auto leaf : LeafDictIdx) {
+                kerr->InsertTuple1( leaf.first, leaf.second );
+            }
+
+            auto kri = oct->GetLeafData()->AddArray(kr);
+            auto kii = oct->GetLeafData()->AddArray(ki);
+            auto kmi = oct->GetLeafData()->AddArray(km);
+            auto kidi = oct->GetLeafData()->AddArray(kid);
+            auto keri = oct->GetLeafData()->AddArray(kerr);
+
+            auto write = vtkXMLHyperOctreeWriter::New();
+                //write.SetDataModeToAscii()
+                write->SetInputData(oct);
+                std::string fname = std::string("octree-") + to_string(ilay)
+                                  + std::string("-") + to_string(iq) + std::string(".vto");
+                write->SetFileName(fname.c_str());
+                write->Write();
+                write->Delete();
+
+            oct->GetLeafData()->RemoveArray( kri );
+            oct->GetLeafData()->RemoveArray( kii );
+            oct->GetLeafData()->RemoveArray( kmi );
+            oct->GetLeafData()->RemoveArray( kidi );
+            oct->GetLeafData()->RemoveArray( keri );
+
+            kerr->Delete();
+            kid->Delete();
+            kr->Delete();
+            ki->Delete();
+            km->Delete();
+
+            }
+
+            curse->Delete();
+            oct->Delete();
+        #else
+            throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
+        #endif
+
+        }
+        std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" <<  Size(0)*Size(1)*Size(2)
+                  << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) <<  std::endl;
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  f
+    //--------------------------------------------------------------------------------------
+    Complex Coupling::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
+        return volume*Ht.dot(Hr);
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  EvaluateKids
+    //--------------------------------------------------------------------------------------
+    void Coupling::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
+        const Complex& parentVal ) {
+
+        std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
+        //std::cout.flush();
+
+        // Next level step, interested in one level below
+        // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
+        Vector3r step  = size.array() / (Real)(1 << (level+1) );
+        Real vol = (step(0)*step(1)*step(2));     // volume of each child
+
+        Vector3r pos =  cpos - step/2.;
+        Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
+                        0,       0,       0,
+                  step[0],       0,       0,
+                        0, step[1],       0,
+                  step[0], step[1],       0,
+                        0,       0, step[2],
+                  step[0],       0, step[2],
+                        0, step[1], step[2],
+                  step[0], step[1], step[2] ).finished();
+
+        VectorXcr kvals(8);       // individual kernel vals
+        cpoints->ClearFields();
+        for (int ichild=0; ichild<8; ++ichild) {
+            Vector3r cp = pos;    // Eigen complains about combining these
+            cp += posadd.row(ichild);
+            cpoints->SetLocation( ichild, cp );
+        }
+
+        Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
+        Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
+        for ( auto EMCalc : EMEarths ) {
+            EMCalc.second->GetFieldPoints()->ClearFields();
+            EMCalc.second->CalculateWireAntennaFields();
+            switch (EMCalc.second->GetTxRxMode()) {
+                case TX:
+                    Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    break;
+                case RX:
+                    Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    break;
+                case TXRX:
+                    Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    break;
+                default:
+                    break;
+            }
+        }
+
+        for (int ichild=0; ichild<8; ++ichild) {
+            Vector3r cp = pos;    // Eigen complains about combining these
+            cp += posadd.row(ichild);
+            kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
+        }
+
+        Complex ksum = kvals.sum();     // Kernel sum
+        // Evaluate whether or not furthur splitting is needed
+        if ( std::abs(ksum-parentVal) > tol || level < minLevel && level < maxLevel ) {
+            // Not a leaf dive further in
+            for (int ichild=0; ichild<8; ++ichild) {
+                Vector3r cp = pos; // Eigen complains about combining these
+                cp += posadd.row(ichild);
+                EvaluateKids( size, level+1, cp, kvals(ichild) );
+            }
+            return; // not leaf
+        }
+        // implicit else, is a leaf
+        SUM += ksum;
+        VOLSUM  += 8.*vol;
+        nleaves += 1; // could say += 8 just as fairly
+        return;     // is leaf
+    }
+
+    #ifdef LEMMAUSEVTK
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
+    //--------------------------------------------------------------------------------------
+    void Coupling::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
+        const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
+
+        std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
+        std::cout.flush();
+
+        // Next level step, interested in one level below
+        // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
+        Vector3r step  = size.array() / (Real)(1 << (level+1) );
+        Real vol = (step(0)*step(1)*step(2));         // volume of each child
+
+        Vector3r pos =  cpos - step/2.;
+        Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
+                        0,       0,       0,
+                  step[0],       0,       0,
+                        0, step[1],       0,
+                  step[0], step[1],       0,
+                        0,       0, step[2],
+                  step[0],       0, step[2],
+                        0, step[1], step[2],
+                  step[0], step[1], step[2] ).finished();
+
+        VectorXcr kvals(8);       // individual kernel vals
+        cpoints->ClearFields();
+        for (int ichild=0; ichild<8; ++ichild) {
+            Vector3r cp = pos;    // Eigen complains about combining these
+            cp += posadd.row(ichild);
+            cpoints->SetLocation( ichild, cp );
+        }
+
+        Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
+        Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
+        for ( auto EMCalc : EMEarths ) {
+            //EMCalc->GetFieldPoints()->ClearFields();
+            EMCalc.second->CalculateWireAntennaFields();
+            switch (EMCalc.second->GetTxRxMode()) {
+                case TX:
+                    Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    break;
+                case RX:
+                    Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    break;
+                case TXRX:
+                    Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
+                    break;
+                default:
+                    break;
+            }
+        }
+
+        for (int ichild=0; ichild<8; ++ichild) {
+            Vector3r cp = pos;    // Eigen complains about combining these
+            cp += posadd.row(ichild);
+            kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
+        }
+
+        Complex ksum = kvals.sum();     // Kernel sum
+        // Evaluate whether or not furthur splitting is needed
+        if ( std::abs(ksum-parentVal) > tol || level < minLevel && level < maxLevel ) {
+            oct->SubdivideLeaf(curse);
+            for (int ichild=0; ichild<8; ++ichild) {
+                curse->ToChild(ichild);
+                Vector3r cp = pos; // Eigen complains about combining these
+                cp += posadd.row(ichild);
+                /* Test for position via alternative means */
+                /*
+                Real p[3];
+                GetPosition(curse, p);
+                if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
+                    std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
+                              << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
+                    throw std::runtime_error("doom");
+                }
+                */
+                /* End of position test */
+                EvaluateKids2( size, level+1, cp, kvals.row(ichild), oct, curse );
+                curse->ToParent();
+            }
+            return;  // not a leaf
+        }
+        LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
+        LeafDictIdx[curse->GetLeafId()] = nleaves;
+        Kern.row(ilay) += ksum;
+        VOLSUM += 8*vol;
+        nleaves += 1;
+        return;     // is a leaf
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  Coupling
+    //      Method:  GetPosition
+    //--------------------------------------------------------------------------------------
+    void Coupling::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
+        Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
+        //step  = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
+        p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
+        p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
+        p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
+    }
+
+    #endif
+
+} // ----  end of namespace Lemma  ----
+
+/* vim: set tabstop=4 expandtab */
+/* vim: set filetype=cpp */
+

src/KernelV0.cpp

@@ -345,7 +345,7 @@ namespace Lemma {
         const VectorXcr& parentVal ) {
 
         std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
-        std::cout.flush();
+        //std::cout.flush();
 
         // Next level step, interested in one level below
         // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
@@ -373,9 +373,8 @@ namespace Lemma {
 
         Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
         Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
-        //Eigen::Matrix< Complex, 8, 3 > Bt;
         for ( auto EMCalc : EMEarths ) {
-            //EMCalc->GetFieldPoints()->ClearFields();
+            EMCalc.second->GetFieldPoints()->ClearFields();
             EMCalc.second->CalculateWireAntennaFields();
             switch (EMCalc.second->GetTxRxMode()) {
                 case TX: