More progress, octree grids are working, as are EM calculations on them. Need to add logic for kernel caluclation for a prism, and then depth layers.

CMakeLists.txt

@@ -11,6 +11,12 @@ set_target_properties(merlin PROPERTIES
 # Linking
 target_link_libraries(merlin "lemmacore" "fdem1d" )
 
+# Linking
+if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT ) 
+	target_link_libraries(merlin ${VTK_LIBRARIES})
+#	target_link_libraries(lemmacore "matplot")
+endif()
+
 # Testing
 if (LEMMA_ENABLE_TESTING)
 	add_subdirectory(testing)

examples/CMakeLists.txt

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

examples/KernelV0.cpp

@@ -30,8 +30,8 @@ int main() {
 		earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
 
     // Transmitter loops
-    auto Tx1 = CircularLoop(60, 15, 0, 0);
-    auto Tx2 = CircularLoop(60, 15, 15, 0);
+    auto Tx1 = CircularLoop(60, 15, 50, 50);
+    auto Tx2 = CircularLoop(60, 15, 55, 50);
     //auto Tx1 = CircularLoop(60, 15, 0, 0); // was 60
 
     auto Kern = KernelV0::NewSP();
@@ -44,7 +44,7 @@ int main() {
     // may be more natural to work with?
     std::vector<std::string> tx = {std::string("Coil 1")};
     std::vector<std::string> rx = {std::string("Coil 1")};
-    Kern->CalculateK0( tx, rx );
+    Kern->CalculateK0( tx, rx , false ); //, false );
     //Kern->CalculateK0( "Coil 1", "Coil 1" );
 
 }

include/KernelV0.h

@@ -23,6 +23,13 @@
 #include "LemmaObject.h"
 #include "LayeredEarthEM.h"
 #include "PolygonalWireAntenna.h"
+#include "EMEarth1D.h"
+
+#ifdef LEMMAUSEVTK
+#include "vtkHyperOctree.h"
+#include "vtkHyperOctreeCursor.h"
+#include "vtkXMLHyperOctreeWriter.h"
+#endif
 
 namespace Lemma {
 
@@ -123,9 +130,15 @@ namespace Lemma {
         /**
          *   Calculates a single imaging kernel, however, phased arrays are supported
          *   so that more than one transmitter and/or receiver can be specified.
-         *   @param[in]
+         *   @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.
          */
-        void CalculateK0 (const std::vector< std::string >& tx, const std::vector< std::string >& rx );
+        void CalculateK0 (const std::vector< std::string >& tx, const std::vector< std::string >& rx,
+                bool vtkOutput=false );
 
         // ====================  INQUIRY       =======================
         /**
@@ -148,9 +161,9 @@ namespace Lemma {
         /**
          *  Returns the kernel value for an input prism
          */
-        Complex f( const Vector3r& r, const Real& volume );
+        Complex f( const Vector3r& r, const Real& volume , const Vector3cr& Bt);
 
-        void IntegrateOnOctreeGrid( const Real& tolerance );
+        void IntegrateOnOctreeGrid( const Real& tolerance , bool vtkOutput=false );
 
         /**
          *  Recursive call to integrate a function on an adaptive Octree Grid.
@@ -165,18 +178,33 @@ namespace Lemma {
         bool 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.
+         */
+        bool EvaluateKids2(  const Vector3r& size, const int& level, const Vector3r& cpos,
+                            const Complex& parentVal, vtkHyperOctree* octree, vtkHyperOctreeCursor* curse );
+        #endif
+
         // ====================  DATA MEMBERS  =========================
 
-        Complex                                  SUM;
+        Complex                                   SUM;
 
-        Real                                     tol=1e-3;
+        Real                                      tol=1e-3;
 
-        int                                      nleaves;
+        int                                       nleaves;
 
-        std::shared_ptr< LayeredEarthEM >        SigmaModel = nullptr;
+        std::shared_ptr< LayeredEarthEM >         SigmaModel = nullptr;
 
         std::map< std::string , std::shared_ptr< PolygonalWireAntenna > >  TxRx;
 
+        std::vector< std::shared_ptr<EMEarth1D> > EMEarths;
+
+        #ifdef LEMMAUSEVTK
+        std::map< int, Complex  >                 LeafDict;
+        #endif
+
         /** ASCII string representation of the class name */
         static constexpr auto CName = "KernelV0";
 

src/KernelV0.cpp

@@ -19,7 +19,6 @@
 
 
 #include "KernelV0.h"
-#include "EMEarth1D.h"
 #include "FieldPoints.h"
 
 namespace Lemma {
@@ -103,30 +102,23 @@ namespace Lemma {
     //       Class:  KernelV0
     //      Method:  DeSerialize
     //--------------------------------------------------------------------------------------
-    void KernelV0::CalculateK0 (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx ) {
+    void KernelV0::CalculateK0 (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
+            bool vtkOutput ) {
 
+        // All EM calculations will share same field points
+        auto points = FieldPoints::NewSP();
+            points->SetNumberOfPoints(8);
         for (auto tx : Tx) {
             // Set up EMEarth
-            auto EmEarth = EMEarth1D::NewSP();
-                EmEarth->AttachWireAntenna(TxRx[tx]);
-                EmEarth->AttachLayeredEarthEM(SigmaModel);
-         		EmEarth->SetFieldsToCalculate(H);
+            EMEarths.push_back( EMEarth1D::NewSP() );
+                EMEarths.back()->AttachWireAntenna(TxRx[tx]);
+                EMEarths.back()->AttachLayeredEarthEM(SigmaModel);
+                EMEarths.back()->AttachFieldPoints( points );
+         		EMEarths.back()->SetFieldsToCalculate(H);
                 // TODO query for method, altough with flat antennae, this is fastest
-                EmEarth->SetHankelTransformMethod(ANDERSON801);
-
-        IntegrateOnOctreeGrid( 1e-2 );
-
-// 		EmEarth->AttachFieldPoints(receivers);
-//         //EmEarth->SetHankelTransformMethod(FHTKEY101);
-// 	    EmEarth->CalculateWireAntennaFields();
-//         Vector3Xcr Rx1 = receivers->GetHfield(0);
-//         //receivers->ClearFields();
-//
-// 		//EmEarth->AttachWireAntenna(Tx2);
-// 	    EmEarth->CalculateWireAntennaFields();
-//         Rx1 += receivers->GetHfield(0);
-
+                EMEarths.back()->SetHankelTransformMethod(ANDERSON801);
         }
+        IntegrateOnOctreeGrid( 1e-7, vtkOutput );
 
     }
 
@@ -134,19 +126,43 @@ namespace Lemma {
     //       Class:  KernelV0
     //      Method:  IntegrateOnOctreeGrid
     //--------------------------------------------------------------------------------------
-    void KernelV0::IntegrateOnOctreeGrid( const Real& tolerance) {
+    void KernelV0::IntegrateOnOctreeGrid( const Real& tolerance, bool vtkOutput) {
 
         this->tol = tolerance;
         Vector3r                Size;
             Size << 100,100,100;
-        //Vector3r                Origin;
+        Vector3r                Origin;
+            Origin << 0,0,0;
         Vector3r                cpos;
             cpos << 50,50,50;
         int                     maxlevel;
 
         SUM = 0;
         nleaves = 0;
-        EvaluateKids( Size, 0, cpos, -1e2 );
+        if (!vtkOutput) {
+            EvaluateKids( Size, 0, cpos, 1e6 );
+        } 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, 1e6, oct, curse );
+            auto write = vtkXMLHyperOctreeWriter::New();
+                //write.SetDataModeToAscii()
+                write->SetInputData(oct);
+                write->SetFileName("octree.vto");
+                write->Write();
+                write->Delete();
+            curse->Delete();
+            oct->Delete();
+        #else
+            throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
+        #endif
+
+        }
         std::cout << "SUM\t" << SUM << "\t" << 100*100*100 << "\t" << SUM - Complex(100.*100.*100.) <<  std::endl;
         std::cout << "nleaves\t" << nleaves << std::endl;
 
@@ -156,8 +172,10 @@ namespace Lemma {
     //       Class:  KernelV0
     //      Method:  f
     //--------------------------------------------------------------------------------------
-    Complex KernelV0::f( const Vector3r& r, const Real& volume ) {
-        return Complex(volume);
+    Complex KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Bt ) {
+        //std::cout << volume*Bt.norm() << std::endl;
+        return Complex(volume*Bt.norm());
+        //return Complex(volume);
     }
 
     //--------------------------------------------------------------------------------------
@@ -184,31 +202,118 @@ namespace Lemma {
                         0, step[1], step[2],
                   step[0], step[1], step[2] ).finished();
 
+        VectorXcr kvals(8);       // individual kernel vals
+        FieldPoints* cpoints = EMEarths[0]->GetFieldPoints();
+            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 );
+        }
+
+        Vector3Xcr Bt;
+        //Eigen::Matrix< Complex, 8, 3 > Bt;
+        for ( auto EMCalc : EMEarths ) {
+            //EMCalc->GetFieldPoints()->ClearFields();
+            EMCalc->CalculateWireAntennaFields();
+            Bt = EMCalc->GetFieldPoints()->GetHfield(0);
+        }
+
+        for (int ichild=0; ichild<8; ++ichild) {
+            Vector3r cp = pos;    // Eigen complains about combining these
+            cp += posadd.row(ichild);
+            kvals(ichild) = f(cp, vol, Bt.col(ichild));
+        }
+
+        Complex ksum = kvals.sum();     // Kernel sum
+        // Evaluate whether or not furthur splitting is needed
+        if ( std::abs(ksum - parentVal) > tol || level < 5 ) {
+            for (int ichild=0; ichild<8; ++ichild) {
+                Vector3r cp = pos; // Eigen complains about combining these
+                cp += posadd.row(ichild);
+                bool isleaf = EvaluateKids( size, level+1, cp, kvals(ichild) );
+                if (isleaf) {      // Include result in final integral
+                    SUM += ksum;
+                    nleaves += 1;
+                }
+            }
+            return false;  // not leaf
+        }
+        // Save here instead?
+        return true;       // leaf
+    }
+
+    #ifdef LEMMAUSEVTK
+    //--------------------------------------------------------------------------------------
+    //       Class:  KernelV0
+    //      Method:  EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
+    //--------------------------------------------------------------------------------------
+    bool KernelV0::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
+        const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
+
+        std::cout << "\rlevel " << level << "\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+2) );
+
+        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
+        FieldPoints* cpoints = EMEarths[0]->GetFieldPoints();
+            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 );
+        }
+
+        Vector3Xcr Bt;
+        for ( auto EMCalc : EMEarths ) {
+            //EMCalc->GetFieldPoints()->ClearFields();
+            EMCalc->CalculateWireAntennaFields();
+            Bt = EMCalc->GetFieldPoints()->GetHfield(0);
+        }
+
         for (int ichild=0; ichild<8; ++ichild) {
             Vector3r cp = pos; // Eigen complains about combining these
             cp += posadd.row(ichild);
-            kvals(ichild) = f(cp, vol);
+            kvals(ichild) = f(cp, vol, Bt.col(ichild));
         }
-        Complex ksum = kvals.sum();     // Kernel sum
 
+        Complex ksum = kvals.sum();     // Kernel sum
         // Evaluate whether or not furthur splitting is needed
-        if ( std::abs(ksum - parentVal) > tol || level < 5 ) {
+        if ( std::abs(ksum - parentVal) > tol || level < 3 ) {
+            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);
-                bool isleaf = EvaluateKids( size, level+1, pos, kvals(ichild) );
+                bool isleaf = EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
                 if (isleaf) {  // Include result in final integral
-//                  Id = curse.GetLeafId()     // VTK
-//                  LeafDict[Id] = vals[child] // VTK
+                    LeafDict[curse->GetLeafId()] = kvals(ichild);       // VTK
                     SUM += ksum;
                     nleaves += 1;
                 }
+                curse->ToParent();
             }
             return false;  // not leaf
         }
         return true;       // leaf
     }
+    #endif
 
 } // ----  end of namespace Lemma  ----