Few changes to coupling code

examples/Coupling.cpp

@@ -44,13 +44,13 @@ int main(int argc, char** argv) {
         Kern->PushCoil( "Coil 2", Tx2 );
         Kern->SetLayeredEarthEM( earth );
 
-        Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
-        Kern->SetIntegrationOrigin( (Vector3r() << 0,0,2).finished() );
+        Kern->SetIntegrationSize( (Vector3r() << 200,200,50).finished() );
+        Kern->SetIntegrationOrigin( (Vector3r() << 0,0,5).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
+    VectorXr Offsets = VectorXr::LinSpaced(60, 5.25, 60); // nbins, low, high
 
     auto outfile = std::ofstream("coupling.dat");
     for (int ii=0; ii< Offsets.size(); ++ii) {

include/Coupling.h

@@ -213,7 +213,7 @@ namespace Lemma {
          *  HyperOctree, which is useful for visualization.
          */
         void EvaluateKids2(  const Vector3r& size, const int& level, const Vector3r& cpos,
-                            const VectorXcr& parentVal, vtkHyperOctree* octree, vtkHyperOctreeCursor* curse );
+                            const Complex& parentVal, vtkHyperOctree* octree, vtkHyperOctreeCursor* curse );
 
         void GetPosition( vtkHyperOctreeCursor* Cursor, Real* p );
         #endif
@@ -241,7 +241,7 @@ namespace Lemma {
         std::map< std::string , std::shared_ptr< EMEarth1D > >             EMEarths;
 
         #ifdef LEMMAUSEVTK
-        std::map< int, VectorXcr  >               LeafDict;
+        std::map< int, Complex  >                 LeafDict;
         std::map< int, int     >                  LeafDictIdx;
         std::map< int, Real     >                 LeafDictErr;
         #endif

src/Coupling.cpp

@@ -163,9 +163,7 @@ namespace Lemma {
                 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) {
+            EvaluateKids2( Size, 0, cpos, Complex(100.0), oct, curse );
 
             // Fill in leaf data
             vtkDoubleArray* kr = vtkDoubleArray::New();
@@ -190,9 +188,9 @@ namespace Lemma {
 
             //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)) );
+                kr->InsertTuple1( leaf.first, std::real(leaf.second) );
+                ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
+                km->InsertTuple1( leaf.first, std::abs(leaf.second) );
                 kid->InsertTuple1( leaf.first, leaf.first );
                 //LeafVol += std::real(leaf.second);
             }
@@ -211,8 +209,7 @@ namespace Lemma {
             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");
+                std::string fname = std::string("octree-couple") + std::string(".vto");
                 write->SetFileName(fname.c_str());
                 write->Write();
                 write->Delete();
@@ -229,8 +226,6 @@ namespace Lemma {
             ki->Delete();
             km->Delete();
 
-            }
-
             curse->Delete();
             oct->Delete();
         #else
@@ -248,6 +243,7 @@ namespace Lemma {
     //--------------------------------------------------------------------------------------
     Complex Coupling::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
         return volume*Ht.dot(Hr);
+        //return Ht.dot(Hr);
     }
 
     //--------------------------------------------------------------------------------------
@@ -258,7 +254,7 @@ namespace Lemma {
         const Complex& 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)
@@ -410,14 +406,14 @@ namespace Lemma {
                 }
                 */
                 /* End of position test */
-                EvaluateKids2( size, level+1, cp, kvals.row(ichild), oct, curse );
+                EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
                 curse->ToParent();
             }
             return;  // not a leaf
         }
         LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
         LeafDictIdx[curse->GetLeafId()] = nleaves;
-        Kern.row(ilay) += ksum;
+        SUM += ksum;
         VOLSUM += 8*vol;
         nleaves += 1;
         return;     // is a leaf