first working version of particle 2d fmm (Liu-style)

euklid · Aug 22, 2015 · 54f3a67 · 54f3a67
1 parent 1b9f795
commit 54f3a67
Show file tree

Hide file tree

Showing 7 changed files with 75 additions and 13 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -42,4 +42,4 @@ set(HEADERS
             fmm_gmres_solver.h
 )
 add_executable(${PROJECT_NAME} ${SOURCES} ${HEADERS})
-target_link_libraries(${PROJECT_NAME} armadillo superlu)
+target_link_libraries(${PROJECT_NAME} armadillo superlu m)
diff --git a/TESTING/test_complex.cpp b/TESTING/test_complex.cpp
@@ -59,7 +59,6 @@ int main()
 
     assert(b == (b*a)/a);
     c = b;
-
     std::cout << c.real << " " << c.img << std::endl;
     std::cout << b.real << " " << b.img << std::endl;
     std::cout << a.real << " " << a.img << std::endl;
@@ -76,5 +75,10 @@ int main()
     a = complex_t::log(complex_t(-1,1));
     std::cout << a.real <<" " << a.img << std::endl;
 
+    //division
+    a = complex_t(1,2);
+    b = complex_t(3,4);
+    ab =a*b;
+    assert(b == ab/a);
 
 }
diff --git a/cell2d.cpp b/cell2d.cpp
@@ -7,6 +7,9 @@
 
 #include "cell2d.h"
 #include <cmath>
+#ifdef DEBUG
+#include <iostream>
+#endif
 
 
 Cell2D::Cell2D( double size, Point const & center)
@@ -77,6 +80,18 @@ std::vector<Cell*> & Cell2D::divide()
         {
             m_children.push_back(new_cells[i]);
         }
+#ifdef DEBUG
+        unsigned int num_new_cell_tgt_el = new_cells_tgt_elements[i].size();
+        unsigned int num_new_cell_src_el = new_cells_src_elements[i].size();
+        for(unsigned int el = 0 ; el< num_new_cell_src_el; el++)
+        {
+            assert(new_cells[i]->contains_point(new_cells_src_elements[i][el]->get_position()));
+        }
+        for(unsigned int el = 0; el < num_new_cell_tgt_el; el ++)
+        {
+            assert(new_cells[i]->contains_point(new_cells_tgt_elements[i][el]->get_position()));
+        }
+#endif
     }
     return m_children;
     /*for (int i = 1; i>-2; i=i-2)
@@ -143,7 +158,7 @@ bool Cell2D::contains_point(Point const &pt) const
     {
         double abs_diff = diff[i];
         abs_diff = (abs_diff>0)?abs_diff:-abs_diff;
-        if (std::abs(diff[i]) >  half_size)
+        if (abs_diff >  half_size)
         {
             return false;
         }

diff --git a/fmm2d.cpp b/fmm2d.cpp
@@ -118,7 +118,7 @@ void add_moments(std::vector<T> const & summand, std::vector<T> & moments)
 
 void FMM2D::m2l_downward_pass(Cell* cur_cell)
 {
-    std::vector<complex_t> shifted_moments, local_exps = cur_cell->get_local_exps_cmp();
+    std::vector<complex_t> local_exps = cur_cell->get_local_exps_cmp();
     if(!local_exps.size())
     {
         local_exps.resize(m_loc_terms,0);
@@ -127,6 +127,7 @@ void FMM2D::m2l_downward_pass(Cell* cur_cell)
     unsigned int num_interaction_list = interaction_list.size();
     for(unsigned int i = 0; i<num_interaction_list; i++)
     {
+        std::vector<complex_t> shifted_moments(m_loc_terms,0);
         m_kernel->M2L_cmp(interaction_list[i]->get_moments_cmp(),
                           interaction_list[i]->get_center(),
                           shifted_moments,
@@ -240,7 +241,7 @@ void FMM2D::downward_pass()
     Tree_Iterator* it = m_tree->downward_iterator();
     // level 2 only M2L and direct, no L2L
     Cell* cur_cell;
-
+    
     if(m_make_prec)
     {
         if(!m_precond.size())
@@ -261,14 +262,49 @@ void FMM2D::downward_pass()
     {
         cur_cell = it->next();
         if(cur_cell->get_level() > 2) break;
+        if(cur_cell->get_target_elements().empty())
+        {
+            continue;
+        }
+#ifdef DEBUG
+    std::cout << "Downward pass: " << cur_cell->debug_info() << std::endl;
+#endif
+
         //M2L
         m2l_downward_pass(cur_cell);
         direct_downward_pass(cur_cell);
+        // leaf at lvl 2
+        if(cur_cell->is_leaf())
+        {
+            evaluate_far_interactions(cur_cell);
+        }
     }
 
+    // there are no level 3 cells
+    if(cur_cell->get_level() < 3) 
+    {
+        return;
+    }
+
     //reached lvl 3 or higher --> use L2L and M2L
     while(1)
     {
+        if(cur_cell->get_target_elements().empty())
+        {
+            // skip this and go to next if available
+            if(it->has_next())
+            {
+                cur_cell = it->next();
+                continue;
+            }
+            else //last element, leave loop
+            {
+                break;
+            }
+        }
+#ifdef DEBUG
+    std::cout << "Downward pass: " << cur_cell->debug_info() << std::endl;
+#endif
         m2l_downward_pass(cur_cell);
         l2l_downward_pass(cur_cell);
 
@@ -282,7 +318,8 @@ void FMM2D::downward_pass()
         if(it->has_next())
         {
             cur_cell = it->next();
-        } else //last element, leave loop
+        } 
+        else //last element, leave loop
         {
             break;
         }

diff --git a/kernel_laplace_point_2d.cpp b/kernel_laplace_point_2d.cpp
@@ -1,4 +1,7 @@
 #include "kernel_laplace_point_2d.h"
+#ifdef DEBUG
+#include <iostream>
+#endif
 
 double Laplace2DKernel::direct(Element const & target, Element const & src) const
 {
@@ -41,6 +44,9 @@ std::vector<complex_t> Laplace2DKernel::calc_moments_cmp(const std::vector<Eleme
         for(int j = 0; j<num_moments; j++)
         {
             moments[j]+= contrib;
+#ifdef DEBUG
+            std::cout << "calc_moments_cmp: moment contrib "<<j<<" " << contrib.real << " " << contrib.img << std::endl;
+#endif
             contrib*=fac/(j+1);
         }
     }
@@ -146,7 +152,6 @@ void Laplace2DKernel::M2L_cmp(std::vector<complex_t> const & moments,
                               std::vector<complex_t> & loc_exp,
                               complex_t const & loc_center) const
 {
-    loc_exp.resize(moments.size(),0);
     // compute (k-1)!/(z_l-z_c)^k for k>0 and -log(z_l-z_c) for k = 0 for reuse
     assert(mom_center != loc_center);
 
@@ -159,7 +164,7 @@ void Laplace2DKernel::M2L_cmp(std::vector<complex_t> const & moments,
         factors[i] = (factors[i-1]*(i-1))/diff;
     }
 
-    int sign = 1;
+    double sign = 1;
     for(int i = 0; i<loc_exp.size(); i++)
     {
         for(int j = 0; j<moments.size(); j++)

diff --git a/main.cpp b/main.cpp
@@ -188,22 +188,23 @@ int main(int argc, char** argv)
     std::cout << "FMM with " << src_elements.size() << " sources and " 
             << tgt_elements.size() << " targets took " << fmm_sec << " seconds and "
             << fmm_n_sec << " nanoseconds" << std::endl;
-#define OUTPUT 1
-#if OUTPUT
+#define OUTPUT_FMM 0
+#if OUTPUT_FMM
 
     for(std::vector<Element*>::const_iterator it = tgt_elements.begin(); 
             it !=tgt_elements.end(); 
             ++it)
     {
         std::cout << (*it)->get_target_value() << std::endl;
     }
-    
+#endif
 
     // direct method to compare
     std::vector<double> direct_val = direct_method(argv[1]);
 
     //calculate errors and print them
-
+#define OUTPUT_COMP 1
+#if OUTPUT_COMP
     std::cout << direct_val.size() << std::endl;
     unsigned int num_tgts = direct_val.size();
     for(int i = 0; i<num_tgts; i++)

diff --git a/tree2d.cpp b/tree2d.cpp
@@ -72,7 +72,7 @@ void Tree2D::generate_cells(int max_elements)
         assert(m_cells.size()-1 == index);
         index++;
 
-	if(cur_cell->number_of_elements() > max_elements)
+	if(cur_cell->number_of_elements() > max_elements || cur_cell->get_level()<2)
         {
             std::vector<Cell*> & new_cells = cur_cell->divide();
             for (int i = 0; i<new_cells.size(); i++)