Cleaning up kernel calculatio

examples/KernelV0.cpp

@@ -49,6 +49,14 @@ int main() {
         //Kern->SetTolerance( .55 ) ; // 1%
 
         Kern->SetPulseDuration(0.020);
+        VectorXr I(36);
+        I << 397.4208916184016, 352.364477036168, 313.0112765842783, 278.37896394065376, 247.81424224324982,
+             220.77925043190442, 196.76493264105017, 175.31662279234038, 156.0044839325404, 138.73983004230124,
+             123.42064612625474, 109.82713394836259, 97.76534468972267, 87.06061858367781, 77.56000002944572, 69.1280780096311,
+             61.64250263640252, 54.99473044877554, 49.091182970515476, 43.84634004556388, 39.184136917167976, 35.03619319797924,
+             31.347205894128976, 28.06346770557137, 25.139117042424758, 22.53420773366429, 20.214205433283347,
+             18.144318026099942, 16.299965972298878, 14.652633628829891, 13.184271405688083, 11.870540177313893,
+             10.697057141915716, 9.64778948429609, 8.709338689612677, 7.871268012862094;
         Kern->SetPulseCurrent( VectorXr::LinSpaced( 1, 10, 200 )  ); // nbins, low, high
         Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 2, 5, 5.5 ) );
 
@@ -56,7 +64,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, true );
+    Kern->CalculateK0( tx, rx, false );
 
 }
 

include/KernelV0.h

@@ -248,7 +248,7 @@ namespace Lemma {
          *  @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
          */
-        bool EvaluateKids(  const Vector3r& size, const int& level, const Vector3r& cpos,
+        void EvaluateKids(  const Vector3r& size, const int& level, const Vector3r& cpos,
                             const Complex& parentVal );
 
         #ifdef LEMMAUSEVTK
@@ -256,7 +256,7 @@ namespace Lemma {
          *  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,
+        void EvaluateKids2(  const Vector3r& size, const int& level, const Vector3r& cpos,
                             const Complex& parentVal, vtkHyperOctree* octree, vtkHyperOctreeCursor* curse );
 
         void GetPosition( vtkHyperOctreeCursor* Cursor, Real* p );

src/KernelV0.cpp

@@ -337,7 +337,7 @@ namespace Lemma {
     //       Class:  KernelV0
     //      Method:  EvaluateKids
     //--------------------------------------------------------------------------------------
-    bool KernelV0::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
+    void KernelV0::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;
@@ -346,9 +346,7 @@ namespace Lemma {
         // 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) );
-        Vector3r step2 = size.array() / (Real)(1 << (level+2) );
-
-        Real vol = (step2(0)*step2(1)*step2(2));     // volume of each child
+        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>() <<
@@ -403,17 +401,15 @@ namespace Lemma {
             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;
-                    VOLSUM += 8.*vol;
-                    nleaves += 1;
-                }
+                EvaluateKids( size, level+1, cp, kvals(ichild) );
             }
-            return false;  // not leaf
+            return; // not leaf
         }
         // Save here instead?
-        return true;       // leaf
+        SUM += ksum;
+        VOLSUM += 8.*vol;
+        nleaves += 1;
+        return;     // is leaf
     }
 
     #ifdef LEMMAUSEVTK
@@ -421,7 +417,7 @@ namespace Lemma {
     //       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,
+    void KernelV0::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;
@@ -432,9 +428,6 @@ namespace Lemma {
         Vector3r step  = size.array() / (Real)(1 << (level+1) );
         Real vol = (step(0)*step(1)*step(2));         // volume of each child
 
-        Vector3r step2 = size.array() / (Real)(1 << (level+2) );
-        Real vol2 = (step2(0)*step2(1)*step2(2));     // volume of each child
-
         Vector3r pos =  cpos - step/2.;
         Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
                         0,       0,       0,
@@ -483,9 +476,8 @@ namespace Lemma {
 
         Complex ksum = kvals.sum();     // Kernel sum
         // Evaluate whether or not furthur splitting is needed
-        if ( std::abs(ksum-parentVal) > tol ) { // || level < minLevel && level < maxLevel ) {
+        if ( std::abs(ksum - parentVal) > tol || level < minLevel && level < maxLevel ) {
         //if ( std::abs(ksum.real()-parentVal.real())>tol || std::abs(ksum.imag()-parentVal.imag()) >tol ) {
-        //if ( std::abs(ksum) > tol && level < maxLevel ) {
             oct->SubdivideLeaf(curse);
             for (int ichild=0; ichild<8; ++ichild) {
                 curse->ToChild(ichild);
@@ -502,26 +494,17 @@ namespace Lemma {
                 }
                 */
                 /* End of position test */
-                bool isleaf = EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
-                /*
-                if (isleaf) {  // Include result in final integral
-                    LeafDict[curse->GetLeafId()] = ksum/(8.*vol);     //kvals(ichild) / vol;     // VTK
-                    LeafDictIdx[curse->GetLeafId()] = nleaves;        // VTK
-                    SUM += ksum;
-                    VOLSUM += 8*vol;
-                    nleaves += 1;
-                }
-                */
+                EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
                 curse->ToParent();
             }
-            return false;  // not leaf
+            return;  // not a leaf
         }
         LeafDict[curse->GetLeafId()] = ksum/(8.*vol);     //kvals(ichild) / vol;     // VTK
         LeafDictIdx[curse->GetLeafId()] = nleaves;        // VTK
         SUM += ksum;
         VOLSUM += 8*vol;
         nleaves += 1;
-        return true;       // leaf
+        return;     // is a leaf
     }
 
     //--------------------------------------------------------------------------------------