/*
* $Revision: 2552 $
*
* last checkin:
* $Author: gutwenger $
* $Date: 2012-07-05 16:45:20 +0200 (Do, 05. Jul 2012) $
***************************************************************/
/** \file
* \brief Implementation of class NMM (New Multipole Method).
*
* \author Stefan Hachul
*
* \par License:
* This file is part of the Open Graph Drawing Framework (OGDF).
*
* \par
* Copyright (C)
* See README.txt in the root directory of the OGDF installation for details.
*
* \par
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* Version 2 or 3 as published by the Free Software Foundation;
* see the file LICENSE.txt included in the packaging of this file
* for details.
*
* \par
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* \par
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*
* \see http://www.gnu.org/copyleft/gpl.html
***************************************************************/
#include "NMM.h"
#include "../../energybased/FMMMLayout.h"
#include "../../basic/Math.h"
#include "../../energybased/numexcept.h"
#include
#define MIN_BOX_LENGTH 1e-300
#ifdef __BORLANDC__
using _STL::log;
#else
using std::log;
#endif
namespace ogdf {
NMM::NMM()
{
//set MIN_NODE_NUMBER and using_NMM
MIN_NODE_NUMBER = 175; using_NMM = true;
//setting predefined parameters
precision(4); particles_in_leaves(25);
tree_construction_way(FMMMLayout::rtcSubtreeBySubtree);
find_sm_cell(FMMMLayout::scfIteratively);
}
void NMM::calculate_repulsive_forces(
const Graph &G,
NodeArray & A,
NodeArray& F_rep)
{
if(using_NMM) //use NewMultipoleMethod
calculate_repulsive_forces_by_NMM(G,A,F_rep);
else //used the exact naive way
calculate_repulsive_forces_by_exact_method(G,A,F_rep);
}
void NMM::calculate_repulsive_forces_by_NMM(
const Graph &G,
NodeArray& A,
NodeArray& F_rep)
{
QuadTreeNM T;
node v;
DPoint nullpoint (0,0);
NodeArray F_direct(G);
NodeArray F_local_exp(G);
NodeArray F_multipole_exp(G);
List quad_tree_leaves;
//initializations
forall_nodes(v,G)
F_direct[v]=F_local_exp[v]=F_multipole_exp[v]=nullpoint;
quad_tree_leaves.clear();
if(tree_construction_way() == FMMMLayout::rtcPathByPath)
build_up_red_quad_tree_path_by_path(G,A,T);
else //tree_construction_way == FMMMLayout::rtcSubtreeBySubtree
build_up_red_quad_tree_subtree_by_subtree(G,A,T);
form_multipole_expansions(A,T,quad_tree_leaves);
calculate_local_expansions_and_WSPRLS(A,T.get_root_ptr());
transform_local_exp_to_forces(A,quad_tree_leaves,F_local_exp);
transform_multipole_exp_to_forces(A,quad_tree_leaves,F_multipole_exp);
calculate_neighbourcell_forces(A,quad_tree_leaves,F_direct);
add_rep_forces(G,F_direct,F_multipole_exp,F_local_exp,F_rep);
delete_red_quad_tree_and_count_treenodes(T);
}
inline void NMM::calculate_repulsive_forces_by_exact_method(
const Graph &G,
NodeArray& A,
NodeArray& F_rep)
{
ExactMethod.calculate_exact_repulsive_forces(G,A,F_rep);
}
void NMM::make_initialisations(
const Graph &G,
double bl,
DPoint d_l_c,
int p_i_l,
int p,
int t_c_w,
int f_s_c)
{
if(G.numberOfNodes() >= MIN_NODE_NUMBER) //using_NMM
{
using_NMM = true; //indicate that NMM is used for force calculation
particles_in_leaves(p_i_l);
precision(p);
tree_construction_way(t_c_w);
find_sm_cell(f_s_c);
down_left_corner = d_l_c; //Export this two values from FMMM
boxlength = bl;
init_binko(2* precision());
init_power_of_2_array();
}
else //use exact method
{
using_NMM = false; //indicate that exact method is used for force calculation
ExactMethod.make_initialisations(bl,d_l_c,0);
}
}
void NMM::deallocate_memory()
{
if(using_NMM) {
free_binko();
free_power_of_2_array();
}
}
void NMM::update_boxlength_and_cornercoordinate(double b_l, DPoint d_l_c)
{
if(using_NMM) {
boxlength = b_l;
down_left_corner = d_l_c;
}
else
ExactMethod.update_boxlength_and_cornercoordinate(b_l,d_l_c);
}
inline void NMM::init_power_of_2_array()
{
int p = 1;
max_power_of_2_index = 30;
power_of_2 = new int[max_power_of_2_index+1];
for(int i = 0; i<= max_power_of_2_index; i++) {
power_of_2[i] = p;
p*=2;
}
}
inline void NMM::free_power_of_2_array()
{
delete [] power_of_2;
}
inline int NMM::power_of_two(int i)
{
if(i <= max_power_of_2_index)
return power_of_2[i];
else
return static_cast(pow(2.0,i));
}
inline int NMM::maxboxindex (int level)
{
if ((level < 0 )) {
cout <<"Failure NMM::maxboxindex :wrong level "<& A,
QuadTreeNM& T)
{
List act_leaf_List,new_leaf_List;
List *act_leaf_List_ptr,*new_leaf_List_ptr,*help_ptr;
List act_x_List_copy,act_y_List_copy;
QuadTreeNodeNM *act_node_ptr;
build_up_root_node(G,A,T);
act_leaf_List.clear();
new_leaf_List.clear();
act_leaf_List.pushFront(T.get_root_ptr());
act_leaf_List_ptr = &act_leaf_List;
new_leaf_List_ptr = &new_leaf_List;
while(!act_leaf_List_ptr->empty())
{
while(!act_leaf_List_ptr->empty())
{
act_node_ptr = act_leaf_List_ptr->popFrontRet();
make_copy_and_init_Lists(*(act_node_ptr->get_x_List_ptr()),act_x_List_copy,
*(act_node_ptr->get_y_List_ptr()),act_y_List_copy);
T.set_act_ptr(act_node_ptr);
decompose_subtreenode(T,act_x_List_copy,act_y_List_copy,*new_leaf_List_ptr);
}
help_ptr = act_leaf_List_ptr;
act_leaf_List_ptr = new_leaf_List_ptr;
new_leaf_List_ptr = help_ptr;
}
}
void NMM::make_copy_and_init_Lists(
List& L_x_orig,
List& L_x_copy,
List& L_y_orig,
List& L_y_copy)
{
ListIterator origin_x_item,copy_x_item,origin_y_item,copy_y_item,
new_cross_ref_item;
ParticleInfo P_x_orig,P_y_orig,P_x_copy,P_y_copy;
bool L_x_orig_traversed = false;
bool L_y_orig_traversed = false;
L_x_copy.clear();
L_y_copy.clear();
origin_x_item = L_x_orig.begin();
while(!L_x_orig_traversed)
{
//reset values
P_x_orig = *origin_x_item;
P_x_orig.set_subList_ptr(NULL); //clear subList_ptr
P_x_orig.set_copy_item(NULL); //clear copy_item
P_x_orig.unmark(); //unmark this element
P_x_orig.set_tmp_cross_ref_item(NULL);//clear tmp_cross_ref_item
//update L_x_copy
P_x_copy = P_x_orig;
L_x_copy.pushBack(P_x_copy);
//update L_x_orig
P_x_orig.set_copy_item(L_x_copy.rbegin());
*origin_x_item = P_x_orig;
if(origin_x_item != L_x_orig.rbegin())
origin_x_item = L_x_orig.cyclicSucc(origin_x_item);
else
L_x_orig_traversed = true;
}
origin_y_item = L_y_orig.begin();
while(!L_y_orig_traversed)
{
//reset values
P_y_orig = *origin_y_item;
P_y_orig.set_subList_ptr(NULL); //clear subList_ptr
P_y_orig.set_copy_item(NULL); //clear copy_item
P_y_orig.set_tmp_cross_ref_item(NULL);//clear tmp_cross_ref_item
P_y_orig.unmark(); //unmark this element
//update L_x(y)_copy
P_y_copy = P_y_orig;
new_cross_ref_item = (*P_y_orig.get_cross_ref_item()).get_copy_item();
P_y_copy.set_cross_ref_item(new_cross_ref_item);
L_y_copy.pushBack(P_y_copy);
P_x_copy = *new_cross_ref_item;
P_x_copy.set_cross_ref_item(L_y_copy.rbegin());
*new_cross_ref_item = P_x_copy;
//update L_y_orig
P_y_orig.set_copy_item(L_y_copy.rbegin());
*origin_y_item = P_y_orig;
if(origin_y_item != L_y_orig.rbegin())
origin_y_item = L_y_orig.cyclicSucc(origin_y_item);
else
L_y_orig_traversed = true;
}
}
void NMM::build_up_root_node(
const Graph& G,
NodeArray& A,
QuadTreeNM& T)
{
T.init_tree();
T.get_root_ptr()->set_Sm_level(0);
T.get_root_ptr()->set_Sm_downleftcorner(down_left_corner);
T.get_root_ptr()->set_Sm_boxlength(boxlength);
//allocate space for L_x and L_y List of the root node
T.get_root_ptr()->set_x_List_ptr(OGDF_NEW List);
T.get_root_ptr()->set_y_List_ptr(OGDF_NEW List);
create_sorted_coordinate_Lists(G, A, *(T.get_root_ptr()->get_x_List_ptr()), *(T.get_root_ptr()->get_y_List_ptr()));
}
void NMM::create_sorted_coordinate_Lists(
const Graph& G,
NodeArray & A,
List& L_x,
List& L_y)
{
ParticleInfo P_x,P_y;
ListIterator x_item,y_item;
node v;
//build up L_x,L_y and link the Lists
forall_nodes(v,G)
{
P_x.set_x_y_coord(A[v].get_x());
P_y.set_x_y_coord(A[v].get_y());
P_x.set_vertex(v);
P_y.set_vertex(v);
L_x.pushBack(P_x);
L_y.pushBack(P_y);
P_x.set_cross_ref_item(L_y.rbegin());
P_y.set_cross_ref_item(L_x.rbegin());
*L_x.rbegin() = P_x;
*L_y.rbegin() = P_y;
}
//sort L_x and update the links of L_y
ParticleInfoComparer comp;
L_x.quicksort(comp);//Quicksort L_x
for(x_item = L_x.begin(); x_item.valid();++x_item)
{
y_item = (*x_item).get_cross_ref_item();
P_y = *y_item;
P_y.set_cross_ref_item(x_item);
*y_item = P_y;
}
//sort L_y and update the links of L_x
L_y.quicksort(comp);//Quicksort L_x
for(y_item = L_y.begin(); y_item.valid();++y_item)
{
x_item = (*y_item).get_cross_ref_item();
P_x = *x_item;
P_x.set_cross_ref_item(y_item);
*x_item = P_x;
}
}
void NMM::decompose_subtreenode(
QuadTreeNM& T,
List& act_x_List_copy,
List& act_y_List_copy,
List& new_leaf_List)
{
QuadTreeNodeNM* act_ptr = T.get_act_ptr();
int act_particle_number = act_ptr->get_x_List_ptr()->size();
double x_min,x_max,y_min,y_max;
List *L_x_l_ptr,*L_x_r_ptr,*L_x_lb_ptr,*L_x_rb_ptr,*L_x_lt_ptr,
*L_x_rt_ptr;
List *L_y_l_ptr,*L_y_r_ptr,*L_y_lb_ptr,*L_y_rb_ptr,*L_y_lt_ptr,
*L_y_rt_ptr;
L_x_l_ptr = L_x_r_ptr = L_x_lb_ptr = L_x_lt_ptr = L_x_rb_ptr = L_x_rt_ptr = NULL;
L_y_l_ptr = L_y_r_ptr = L_y_lb_ptr = L_y_lt_ptr = L_y_rb_ptr = L_y_rt_ptr = NULL;
calculate_boundaries_of_act_node(T.get_act_ptr(),x_min,x_max,y_min,y_max);
if(find_sm_cell() == FMMMLayout::scfIteratively)
find_small_cell_iteratively(T.get_act_ptr(),x_min,x_max,y_min,y_max);
else //find_small_cell == FMMMLayout::scfAluru
find_small_cell_iteratively(T.get_act_ptr(),x_min,x_max,y_min,y_max);
if( (act_particle_number > particles_in_leaves()) &&
((x_max-x_min >=MIN_BOX_LENGTH) || (y_max-y_min >= MIN_BOX_LENGTH )))
{//if0
//recursive calls for the half of the quad that contains the most particles
split_in_x_direction(act_ptr,L_x_l_ptr,L_y_l_ptr,
L_x_r_ptr,L_y_r_ptr);
if((L_x_r_ptr == NULL) ||
(L_x_l_ptr != NULL && L_x_l_ptr->size() > L_x_r_ptr->size()))
{//if1 left half contains more particles
split_in_y_direction(act_ptr,L_x_lb_ptr,
L_y_lb_ptr,L_x_lt_ptr,L_y_lt_ptr);
if((L_x_lt_ptr == NULL)||
(L_x_lb_ptr != NULL && L_x_lb_ptr->size() > L_x_lt_ptr->size()))
{//if2
T.create_new_lb_child(L_x_lb_ptr,L_y_lb_ptr);
T.go_to_lb_child();
decompose_subtreenode(T,act_x_List_copy,act_y_List_copy,new_leaf_List);
T.go_to_father();
}//if2
else //L_x_lt_ptr != NULL && L_x_lb_ptr->size() <= L_x_lt_ptr->size()
{//else1
T.create_new_lt_child(L_x_lt_ptr,L_y_lt_ptr);
T.go_to_lt_child();
decompose_subtreenode(T,act_x_List_copy,act_y_List_copy,new_leaf_List);
T.go_to_father();
}//else1
}//if1
else //L_x_r_ptr != NULL && (L_x_l_ptr->size() <= L_x_r_ptr->size())
{//else2 right half contains more particles
split_in_y_direction(act_ptr,L_x_rb_ptr,
L_y_rb_ptr,L_x_rt_ptr,L_y_rt_ptr);
if ((L_x_rt_ptr == NULL) ||
(L_x_rb_ptr != NULL && L_x_rb_ptr->size() > L_x_rt_ptr->size()))
{//if3
T.create_new_rb_child(L_x_rb_ptr,L_y_rb_ptr);
T.go_to_rb_child();
decompose_subtreenode(T,act_x_List_copy,act_y_List_copy,new_leaf_List);
T.go_to_father();
}//if3
else// L_x_rt_ptr != NULL && L_x_rb_ptr->size() <= L_x_rt_ptr->size()
{//else3
T.create_new_rt_child(L_x_rt_ptr,L_y_rt_ptr);
T.go_to_rt_child();
decompose_subtreenode(T,act_x_List_copy,act_y_List_copy,new_leaf_List);
T.go_to_father();
}//else3
}//else2
//build up the rest of the quad-subLists
if( L_x_l_ptr != NULL && L_x_lb_ptr == NULL && L_x_lt_ptr == NULL &&
!act_ptr->child_lb_exists() && !act_ptr->child_lt_exists() )
split_in_y_direction(act_ptr,L_x_l_ptr,L_x_lb_ptr,L_x_lt_ptr,L_y_l_ptr,
L_y_lb_ptr,L_y_lt_ptr);
else if( L_x_r_ptr != NULL && L_x_rb_ptr == NULL && L_x_rt_ptr == NULL &&
!act_ptr->child_rb_exists() && !act_ptr->child_rt_exists() )
split_in_y_direction(act_ptr,L_x_r_ptr,L_x_rb_ptr,L_x_rt_ptr,L_y_r_ptr,
L_y_rb_ptr,L_y_rt_ptr);
//create rest of the childnodes
if((!act_ptr->child_lb_exists()) && (L_x_lb_ptr != NULL))
{
T.create_new_lb_child(L_x_lb_ptr,L_y_lb_ptr);
T.go_to_lb_child();
new_leaf_List.pushBack(T.get_act_ptr());
T.go_to_father();
}
if((!act_ptr->child_lt_exists()) && (L_x_lt_ptr != NULL))
{
T.create_new_lt_child(L_x_lt_ptr,L_y_lt_ptr);
T.go_to_lt_child();
new_leaf_List.pushBack(T.get_act_ptr());
T.go_to_father();
}
if((!act_ptr->child_rb_exists()) && (L_x_rb_ptr != NULL))
{
T.create_new_rb_child(L_x_rb_ptr,L_y_rb_ptr);
T.go_to_rb_child();
new_leaf_List.pushBack(T.get_act_ptr());
T.go_to_father();
}
if((!act_ptr->child_rt_exists()) && (L_x_rt_ptr != NULL))
{
T.create_new_rt_child(L_x_rt_ptr,L_y_rt_ptr);
T.go_to_rt_child();
new_leaf_List.pushBack(T.get_act_ptr());
T.go_to_father();
}
//reset act_ptr->set_x(y)_List_ptr to avoid multiple deleting of dynamic memory;
//(only if *act_ptr is a leaf of T the reserved space is freed (and this is
//sufficient !!!))
act_ptr->set_x_List_ptr(NULL);
act_ptr->set_y_List_ptr(NULL);
}//if0
else
{ //else a leaf or machineprecision is reached:
//The List contained_nodes is set for *act_ptr and the information of
//act_x_List_copy and act_y_List_copy is used to insert particles into the
//shorter Lists of previous touched treenodes;additionaly the dynamical allocated
//space for *act_ptr->get_x(y)_List_ptr() is freed.
List L;
ListIterator it;
//set List contained nodes
L.clear();
for(it = act_ptr->get_x_List_ptr()->begin();it.valid();++it)
L.pushBack((*it).get_vertex());
T.get_act_ptr()->set_contained_nodes(L);
//insert particles into previous touched Lists
build_up_sorted_subLists(act_x_List_copy,act_y_List_copy);
//free allocated space for *act_ptr->get_x(y)_List_ptr()
act_ptr->get_x_List_ptr()->clear();//free used space for old L_x,L_y Lists
act_ptr->get_y_List_ptr()->clear();
}//else
}
inline void NMM::calculate_boundaries_of_act_node(
QuadTreeNodeNM* act_ptr,
double& x_min,
double& x_max,
double& y_min,
double& y_max)
{
List* L_x_ptr = act_ptr->get_x_List_ptr();
List* L_y_ptr = act_ptr->get_y_List_ptr();
x_min = (*L_x_ptr->begin()).get_x_y_coord();
x_max = (*L_x_ptr->rbegin()).get_x_y_coord();
y_min = (*L_y_ptr->begin()).get_x_y_coord();
y_max = (*L_y_ptr->rbegin()).get_x_y_coord();
}
bool NMM::in_lt_quad(
QuadTreeNodeNM* act_ptr,
double x_min,
double x_max,
double y_min,
double y_max)
{
double l = act_ptr->get_Sm_downleftcorner().m_x;
double r = act_ptr->get_Sm_downleftcorner().m_x+act_ptr->get_Sm_boxlength()/2;
double b = act_ptr->get_Sm_downleftcorner().m_y+act_ptr->get_Sm_boxlength()/2;
double t = act_ptr->get_Sm_downleftcorner().m_y+act_ptr->get_Sm_boxlength();
if(l <= x_min && x_max < r && b <= y_min && y_max < t )
return true;
else if(x_min == x_max && y_min == y_max && l==r && t == b && x_min ==r && y_min ==b)
return true;
else
return false;
}
bool NMM::in_rt_quad(
QuadTreeNodeNM* act_ptr,
double x_min,
double x_max,
double y_min,
double y_max)
{
double l = act_ptr->get_Sm_downleftcorner().m_x+act_ptr->get_Sm_boxlength()/2;
double r = act_ptr->get_Sm_downleftcorner().m_x+act_ptr->get_Sm_boxlength();
double b = act_ptr->get_Sm_downleftcorner().m_y+act_ptr->get_Sm_boxlength()/2;
double t = act_ptr->get_Sm_downleftcorner().m_y+act_ptr->get_Sm_boxlength();
if(l <= x_min && x_max < r && b <= y_min && y_max < t )
return true;
else if(x_min == x_max && y_min == y_max && l==r && t == b && x_min ==r && y_min ==b)
return true;
else
return false;
}
bool NMM::in_lb_quad(
QuadTreeNodeNM* act_ptr,
double x_min,
double x_max,
double y_min,
double y_max)
{
double l = act_ptr->get_Sm_downleftcorner().m_x;
double r = act_ptr->get_Sm_downleftcorner().m_x+act_ptr->get_Sm_boxlength()/2;
double b = act_ptr->get_Sm_downleftcorner().m_y;
double t = act_ptr->get_Sm_downleftcorner().m_y+act_ptr->get_Sm_boxlength()/2;
if(l <= x_min && x_max < r && b <= y_min && y_max < t )
return true;
else if(x_min == x_max && y_min == y_max && l==r && t == b && x_min ==r && y_min ==b)
return true;
else
return false;
}
bool NMM::in_rb_quad(
QuadTreeNodeNM* act_ptr,
double x_min,
double x_max,
double y_min,
double y_max)
{
double l = act_ptr->get_Sm_downleftcorner().m_x+act_ptr->get_Sm_boxlength()/2;
double r = act_ptr->get_Sm_downleftcorner().m_x+act_ptr->get_Sm_boxlength();
double b = act_ptr->get_Sm_downleftcorner().m_y;
double t = act_ptr->get_Sm_downleftcorner().m_y+act_ptr->get_Sm_boxlength()/2;
if(l <= x_min && x_max < r && b <= y_min && y_max < t )
return true;
else if(x_min == x_max && y_min == y_max && l==r && t == b && x_min ==r && y_min ==b)
return true;
else
return false;
}
void NMM::split_in_x_direction(
QuadTreeNodeNM* act_ptr,
List *& L_x_left_ptr,
List*& L_y_left_ptr,
List *& L_x_right_ptr,
List*& L_y_right_ptr)
{
ListIterator l_item = act_ptr->get_x_List_ptr()->begin();
ListIterator r_item = act_ptr->get_x_List_ptr()->rbegin();
ListIterator last_left_item;
double act_Sm_boxlength_half = act_ptr->get_Sm_boxlength()/2;
double x_mid_coord = act_ptr->get_Sm_downleftcorner().m_x+ act_Sm_boxlength_half;
double l_xcoord,r_xcoord;
bool last_left_item_found = false;
bool left_particleList_empty = false;
bool right_particleList_empty = false;
bool left_particleList_larger = true;
//traverse *act_ptr->get_x_List_ptr() from left and right
while(!last_left_item_found)
{//while
l_xcoord = (*l_item).get_x_y_coord();
r_xcoord = (*r_item).get_x_y_coord();
if(l_xcoord >= x_mid_coord)
{
left_particleList_larger = false;
last_left_item_found = true;
if(l_item != act_ptr->get_x_List_ptr()->begin())
last_left_item = act_ptr->get_x_List_ptr()->cyclicPred(l_item);
else
left_particleList_empty = true;
}
else if(r_xcoord < x_mid_coord)
{
last_left_item_found = true;
if(r_item != act_ptr->get_x_List_ptr()->rbegin())
last_left_item = r_item;
else
right_particleList_empty = true;
}
if(!last_left_item_found)
{
l_item = act_ptr->get_x_List_ptr()->cyclicSucc(l_item);
r_item = act_ptr->get_x_List_ptr()->cyclicPred(r_item);
}
}//while
//get the L_x(y) Lists of the bigger half (from *act_ptr->get_x(y)_List_ptr))
//and make entries in L_x_copy,L_y_copy for the smaller halfs
if(left_particleList_empty)
{
L_x_left_ptr = NULL;
L_y_left_ptr = NULL;
L_x_right_ptr = act_ptr->get_x_List_ptr();
L_y_right_ptr = act_ptr->get_y_List_ptr();
}
else if(right_particleList_empty)
{
L_x_left_ptr = act_ptr->get_x_List_ptr();
L_y_left_ptr = act_ptr->get_y_List_ptr();
L_x_right_ptr = NULL;
L_y_right_ptr = NULL;
}
else if(left_particleList_larger)
x_delete_right_subLists(act_ptr,L_x_left_ptr,L_y_left_ptr,
L_x_right_ptr,L_y_right_ptr,last_left_item);
else //left particleList is smaller or equal to right particleList
x_delete_left_subLists(act_ptr,L_x_left_ptr,L_y_left_ptr,
L_x_right_ptr,L_y_right_ptr,last_left_item);
}
void NMM::split_in_y_direction(
QuadTreeNodeNM* act_ptr,
List*& L_x_left_ptr,
List*& L_y_left_ptr,
List*& L_x_right_ptr,
List*& L_y_right_ptr)
{
ListIterator l_item = act_ptr->get_y_List_ptr()->begin();
ListIterator r_item = act_ptr->get_y_List_ptr()->rbegin();
ListIterator last_left_item;
double act_Sm_boxlength_half = act_ptr->get_Sm_boxlength()/2;
double y_mid_coord = act_ptr->get_Sm_downleftcorner().m_y+ act_Sm_boxlength_half;
double l_ycoord,r_ycoord;
bool last_left_item_found = false;
bool left_particleList_empty = false;
bool right_particleList_empty = false;
bool left_particleList_larger = true;
//traverse *act_ptr->get_y_List_ptr() from left and right
while(!last_left_item_found)
{//while
l_ycoord = (*l_item).get_x_y_coord();
r_ycoord = (*r_item).get_x_y_coord();
if(l_ycoord >= y_mid_coord)
{
left_particleList_larger = false;
last_left_item_found = true;
if(l_item != act_ptr->get_y_List_ptr()->begin())
last_left_item = act_ptr->get_y_List_ptr()->cyclicPred(l_item);
else
left_particleList_empty = true;
}
else if(r_ycoord < y_mid_coord)
{
last_left_item_found = true;
if(r_item != act_ptr->get_y_List_ptr()->rbegin())
last_left_item = r_item;
else
right_particleList_empty = true;
}
if(!last_left_item_found)
{
l_item = act_ptr->get_y_List_ptr()->cyclicSucc(l_item);
r_item = act_ptr->get_y_List_ptr()->cyclicPred(r_item);
}
}//while
//get the L_x(y) Lists of the bigger half (from *act_ptr->get_x(y)_List_ptr))
//and make entries in L_x_copy,L_y_copy for the smaller halfs
if(left_particleList_empty)
{
L_x_left_ptr = NULL;
L_y_left_ptr = NULL;
L_x_right_ptr = act_ptr->get_x_List_ptr();
L_y_right_ptr = act_ptr->get_y_List_ptr();
}
else if(right_particleList_empty)
{
L_x_left_ptr = act_ptr->get_x_List_ptr();
L_y_left_ptr = act_ptr->get_y_List_ptr();
L_x_right_ptr = NULL;
L_y_right_ptr = NULL;
}
else if(left_particleList_larger)
y_delete_right_subLists(act_ptr,L_x_left_ptr,L_y_left_ptr,
L_x_right_ptr,L_y_right_ptr,last_left_item);
else //left particleList is smaller or equal to right particleList
y_delete_left_subLists(act_ptr,L_x_left_ptr,L_y_left_ptr,
L_x_right_ptr,L_y_right_ptr,last_left_item);
}
void NMM::x_delete_right_subLists(
QuadTreeNodeNM* act_ptr,
List *& L_x_left_ptr,
List*& L_y_left_ptr,
List *& L_x_right_ptr,
List*& L_y_right_ptr,
ListIterator last_left_item)
{
ParticleInfo act_p_info,p_in_L_x_info,p_in_L_y_info,del_p_info;
ListIterator act_item,p_in_L_x_item,p_in_L_y_item,del_item;
bool last_item_reached =false;
L_x_left_ptr = act_ptr->get_x_List_ptr();
L_y_left_ptr = act_ptr->get_y_List_ptr();
L_x_right_ptr = OGDF_NEW List;
L_y_right_ptr = OGDF_NEW List;
act_item = L_x_left_ptr->cyclicSucc(last_left_item);
while(!last_item_reached)
{//while
act_p_info = *act_item;
del_item = act_item;
del_p_info = act_p_info;
//save references for *L_x(y)_right_ptr in L_x(y)_copy
p_in_L_x_item = act_p_info.get_copy_item();
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.set_subList_ptr(L_x_right_ptr);
*p_in_L_x_item = p_in_L_x_info;
p_in_L_y_item = (*act_p_info.get_cross_ref_item()).get_copy_item();
p_in_L_y_info = *p_in_L_y_item;
p_in_L_y_info.set_subList_ptr(L_y_right_ptr);
*p_in_L_y_item = p_in_L_y_info;
if(act_item != L_x_left_ptr->rbegin())
act_item = L_x_left_ptr->cyclicSucc(act_item);
else
last_item_reached = true;
//create *L_x(y)_left_ptr
L_y_left_ptr->del(del_p_info.get_cross_ref_item());
L_x_left_ptr->del(del_item);
}//while
}
void NMM::x_delete_left_subLists(
QuadTreeNodeNM* act_ptr,
List *& L_x_left_ptr,
List*& L_y_left_ptr,
List *& L_x_right_ptr,
List*& L_y_right_ptr,
ListIterator last_left_item)
{
ParticleInfo act_p_info,p_in_L_x_info,p_in_L_y_info,del_p_info;
ListIterator act_item,p_in_L_x_item,p_in_L_y_item,del_item;
bool last_item_reached =false;
L_x_right_ptr = act_ptr->get_x_List_ptr();
L_y_right_ptr = act_ptr->get_y_List_ptr();
L_x_left_ptr = OGDF_NEW List;
L_y_left_ptr = OGDF_NEW List;
act_item = L_x_right_ptr->begin();
while(!last_item_reached)
{//while
act_p_info = *act_item;
del_item = act_item;
del_p_info = act_p_info;
//save references for *L_x(y)_right_ptr in L_x(y)_copy
p_in_L_x_item = act_p_info.get_copy_item();
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.set_subList_ptr(L_x_left_ptr);
*p_in_L_x_item = p_in_L_x_info;
p_in_L_y_item =(*act_p_info.get_cross_ref_item()).get_copy_item();
p_in_L_y_info = *p_in_L_y_item;
p_in_L_y_info.set_subList_ptr(L_y_left_ptr);
*p_in_L_y_item = p_in_L_y_info;
if(act_item != last_left_item)
act_item = L_x_right_ptr->cyclicSucc(act_item);
else
last_item_reached = true;
//create *L_x(y)_right_ptr
L_y_right_ptr->del(del_p_info.get_cross_ref_item());
L_x_right_ptr->del(del_item);
}//while
}
void NMM::y_delete_right_subLists(
QuadTreeNodeNM* act_ptr,
List *& L_x_left_ptr,
List*& L_y_left_ptr,
List *& L_x_right_ptr,
List*& L_y_right_ptr,
ListIterator last_left_item)
{
ParticleInfo act_p_info,p_in_L_x_info,p_in_L_y_info,del_p_info;
ListIterator act_item,p_in_L_x_item,p_in_L_y_item,del_item;
bool last_item_reached =false;
L_x_left_ptr = act_ptr->get_x_List_ptr();
L_y_left_ptr = act_ptr->get_y_List_ptr();
L_x_right_ptr = OGDF_NEW List;
L_y_right_ptr = OGDF_NEW List;
act_item = L_y_left_ptr->cyclicSucc(last_left_item);
while(!last_item_reached)
{//while
act_p_info = *act_item;
del_item = act_item;
del_p_info = act_p_info;
//save references for *L_x(y)_right_ptr in L_x(y)_copy
p_in_L_y_item = act_p_info.get_copy_item();
p_in_L_y_info = *p_in_L_y_item;
p_in_L_y_info.set_subList_ptr(L_y_right_ptr);
*p_in_L_y_item = p_in_L_y_info;
p_in_L_x_item = (*act_p_info.get_cross_ref_item()).get_copy_item();
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.set_subList_ptr(L_x_right_ptr);
*p_in_L_x_item = p_in_L_x_info;
if(act_item != L_y_left_ptr->rbegin())
act_item = L_y_left_ptr->cyclicSucc(act_item);
else
last_item_reached = true;
//create *L_x(y)_left_ptr
L_x_left_ptr->del(del_p_info.get_cross_ref_item());
L_y_left_ptr->del(del_item);
}//while
}
void NMM::y_delete_left_subLists(
QuadTreeNodeNM* act_ptr,
List*& L_x_left_ptr,
List*& L_y_left_ptr,
List *& L_x_right_ptr,
List*& L_y_right_ptr,
ListIterator last_left_item)
{
ParticleInfo act_p_info,p_in_L_x_info,p_in_L_y_info,del_p_info;
ListIterator act_item,p_in_L_x_item,p_in_L_y_item,del_item;
bool last_item_reached =false;
L_x_right_ptr = act_ptr->get_x_List_ptr();
L_y_right_ptr = act_ptr->get_y_List_ptr();
L_x_left_ptr = OGDF_NEW List;
L_y_left_ptr = OGDF_NEW List;
act_item = L_y_right_ptr->begin();
while(!last_item_reached)
{//while
act_p_info = *act_item;
del_item = act_item;
del_p_info = act_p_info;
//save references for *L_x(y)_right_ptr in L_x(y)_copy
p_in_L_y_item = act_p_info.get_copy_item();
p_in_L_y_info = *p_in_L_y_item;
p_in_L_y_info.set_subList_ptr(L_y_left_ptr);
*p_in_L_y_item = p_in_L_y_info;
p_in_L_x_item = (*act_p_info.get_cross_ref_item()).get_copy_item();
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.set_subList_ptr(L_x_left_ptr);
*p_in_L_x_item = p_in_L_x_info;
if(act_item != last_left_item)
act_item = L_y_right_ptr->cyclicSucc(act_item);
else
last_item_reached = true;
//create *L_x(y)_right_ptr
L_x_right_ptr->del(del_p_info.get_cross_ref_item());
L_y_right_ptr->del(del_item);
}//while
}
void NMM::split_in_y_direction(
QuadTreeNodeNM* act_ptr,
List*& L_x_ptr,
List*& L_x_b_ptr,
List*& L_x_t_ptr,
List*& L_y_ptr,
List*& L_y_b_ptr,
List*& L_y_t_ptr)
{
ListIterator l_item = L_y_ptr->begin();
ListIterator r_item = L_y_ptr->rbegin();
ListIterator last_left_item;
double act_Sm_boxlength_half = act_ptr->get_Sm_boxlength()/2;
double y_mid_coord = act_ptr->get_Sm_downleftcorner().m_y+ act_Sm_boxlength_half;
double l_ycoord,r_ycoord;
bool last_left_item_found = false;
bool left_particleList_empty = false;
bool right_particleList_empty = false;
bool left_particleList_larger = true;
//traverse *L_y_ptr from left and right
while(!last_left_item_found)
{//while
l_ycoord = (*l_item).get_x_y_coord();
r_ycoord = (*r_item).get_x_y_coord();
if(l_ycoord >= y_mid_coord)
{
left_particleList_larger = false;
last_left_item_found = true;
if(l_item != L_y_ptr->begin())
last_left_item = L_y_ptr->cyclicPred(l_item);
else
left_particleList_empty = true;
}
else if(r_ycoord < y_mid_coord)
{
last_left_item_found = true;
if(r_item != L_y_ptr->rbegin())
last_left_item = r_item;
else
right_particleList_empty = true;
}
if(!last_left_item_found)
{
l_item = L_y_ptr->cyclicSucc(l_item);
r_item = L_y_ptr->cyclicPred(r_item);
}
}//while
//create *L_x_l(b)_ptr
if(left_particleList_empty)
{
L_x_b_ptr = NULL;
L_y_b_ptr = NULL;
L_x_t_ptr = L_x_ptr;
L_y_t_ptr = L_y_ptr;
}
else if(right_particleList_empty)
{
L_x_b_ptr = L_x_ptr;
L_y_b_ptr = L_y_ptr;
L_x_t_ptr = NULL;
L_y_t_ptr = NULL;
}
else if(left_particleList_larger)
y_move_right_subLists(L_x_ptr,L_x_b_ptr,L_x_t_ptr,L_y_ptr,L_y_b_ptr,L_y_t_ptr,
last_left_item);
else //left particleList is smaller or equal to right particleList
y_move_left_subLists(L_x_ptr,L_x_b_ptr,L_x_t_ptr,L_y_ptr,L_y_b_ptr,L_y_t_ptr,
last_left_item);
}
void NMM::y_move_left_subLists(
List*& L_x_ptr,
List *& L_x_l_ptr,
List*& L_x_r_ptr,
List*& L_y_ptr,
List *& L_y_l_ptr,
List*& L_y_r_ptr,
ListIterator last_left_item)
{
ParticleInfo p_in_L_x_info,p_in_L_y_info;
ListIterator p_in_L_x_item,p_in_L_y_item,del_item;
bool last_item_reached =false;
L_x_r_ptr = L_x_ptr;
L_y_r_ptr = L_y_ptr;
L_x_l_ptr = OGDF_NEW List;
L_y_l_ptr = OGDF_NEW List;
p_in_L_y_item = L_y_r_ptr->begin();
//build up the L_y_Lists and update crossreferences in *L_x_l_ptr
while(!last_item_reached)
{//while
p_in_L_y_info = *p_in_L_y_item;
del_item = p_in_L_y_item;
//create *L_x(y)_l_ptr
L_y_l_ptr->pushBack(p_in_L_y_info);
p_in_L_x_item = p_in_L_y_info.get_cross_ref_item();
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.set_cross_ref_item(L_y_l_ptr->rbegin());
p_in_L_x_info.mark(); //mark this element of the List
*p_in_L_x_item = p_in_L_x_info;
if(p_in_L_y_item != last_left_item)
p_in_L_y_item = L_y_r_ptr->cyclicSucc(p_in_L_y_item);
else
last_item_reached = true;
//create *L_y_r_ptr
L_y_r_ptr->del(del_item);
}//while
//build up the L_x Lists and update crossreferences in *L_y_l_ptr
last_item_reached = false;
p_in_L_x_item = L_x_r_ptr->begin();
while(!last_item_reached)
{//while
del_item = p_in_L_x_item;
if((*del_item).is_marked())
{
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.unmark();
L_x_l_ptr->pushBack(p_in_L_x_info);
p_in_L_y_item = p_in_L_x_info.get_cross_ref_item();
p_in_L_y_info = *p_in_L_y_item;
p_in_L_y_info.set_cross_ref_item(L_x_l_ptr->rbegin());
*p_in_L_y_item = p_in_L_y_info;
}
if(p_in_L_x_item != L_x_r_ptr->rbegin())
p_in_L_x_item = L_x_r_ptr->cyclicSucc(p_in_L_x_item);
else
last_item_reached = true;
//create *L_x_r_ptr
if((*del_item).is_marked())
L_x_r_ptr->del(del_item);
}//while
}
void NMM::y_move_right_subLists(
List*& L_x_ptr,
List *& L_x_l_ptr,
List*& L_x_r_ptr,
List*& L_y_ptr,
List *& L_y_l_ptr,
List*& L_y_r_ptr,
ListIterator last_left_item)
{
ParticleInfo p_in_L_x_info,p_in_L_y_info;
ListIterator p_in_L_x_item,p_in_L_y_item,del_item;
bool last_item_reached =false;
L_x_l_ptr = L_x_ptr;
L_y_l_ptr = L_y_ptr;
L_x_r_ptr = OGDF_NEW List;
L_y_r_ptr = OGDF_NEW List;
p_in_L_y_item = L_y_l_ptr->cyclicSucc(last_left_item);
//build up the L_y_Lists and update crossreferences in *L_x_r_ptr
while(!last_item_reached)
{//while
p_in_L_y_info = *p_in_L_y_item;
del_item = p_in_L_y_item;
//create *L_x(y)_r_ptr
L_y_r_ptr->pushBack(p_in_L_y_info);
p_in_L_x_item = p_in_L_y_info.get_cross_ref_item();
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.set_cross_ref_item(L_y_r_ptr->rbegin());
p_in_L_x_info.mark(); //mark this element of the List
*p_in_L_x_item = p_in_L_x_info;
if(p_in_L_y_item != L_y_l_ptr->rbegin())
p_in_L_y_item = L_y_l_ptr->cyclicSucc(p_in_L_y_item);
else
last_item_reached = true;
//create *L_y_l_ptr
L_y_l_ptr->del(del_item);
}//while
//build up the L_x Lists and update crossreferences in *L_y_r_ptr
last_item_reached = false;
p_in_L_x_item = L_x_l_ptr->begin();
while(!last_item_reached)
{//while
del_item = p_in_L_x_item;
if((*del_item).is_marked())
{
p_in_L_x_info = *p_in_L_x_item;
p_in_L_x_info.unmark();
L_x_r_ptr->pushBack(p_in_L_x_info);
p_in_L_y_item = p_in_L_x_info.get_cross_ref_item();
p_in_L_y_info = *p_in_L_y_item;
p_in_L_y_info.set_cross_ref_item(L_x_r_ptr->rbegin());
*p_in_L_y_item = p_in_L_y_info;
}
if(p_in_L_x_item != L_x_l_ptr->rbegin())
p_in_L_x_item = L_x_l_ptr->cyclicSucc(p_in_L_x_item);
else
last_item_reached = true;
//create *L_x_r_ptr
if((*del_item).is_marked())
L_x_l_ptr->del(del_item);
}//while
}
void NMM::build_up_sorted_subLists(
List& L_x_copy,
List& L_y_copy)
{
ParticleInfo P_x,P_y;
List *L_x_ptr,*L_y_ptr;
ListIterator it,new_cross_ref_item;
for(it = L_x_copy.begin();it.valid();++it)
if((*it).get_subList_ptr() != NULL)
{
//reset values
P_x = *it;
L_x_ptr = P_x.get_subList_ptr();
P_x.set_subList_ptr(NULL); //clear subList_ptr
P_x.set_copy_item(NULL); //clear copy_item
P_x.unmark(); //unmark this element
P_x.set_tmp_cross_ref_item(NULL);//clear tmp_cross_ref_item
//update *L_x_ptr
L_x_ptr->pushBack(P_x);
//update L_x_copy
P_x.set_tmp_cross_ref_item(L_x_ptr->rbegin());
*it = P_x;
}
for(it = L_y_copy.begin();it.valid();++it)
if((*it).get_subList_ptr() != NULL)
{
//reset values
P_y = *it;
L_y_ptr = P_y.get_subList_ptr();
P_y.set_subList_ptr(NULL); //clear subList_ptr
P_y.set_copy_item(NULL); //clear copy_item
P_y.unmark(); //unmark this element
P_y.set_tmp_cross_ref_item(NULL);//clear tmp_cross_ref_item
//update *L_x(y)_ptr
new_cross_ref_item = (*P_y.get_cross_ref_item()).get_tmp_cross_ref_item();
P_y.set_cross_ref_item(new_cross_ref_item);
L_y_ptr->pushBack(P_y);
P_x = *new_cross_ref_item;
P_x.set_cross_ref_item(L_y_ptr->rbegin());
*new_cross_ref_item = P_x;
}
}
// **********Functions needed for subtree by subtree tree construction(Begin)*********
void NMM::build_up_red_quad_tree_subtree_by_subtree(
const Graph& G,
NodeArray& A,
QuadTreeNM& T)
{
List act_subtree_root_List,new_subtree_root_List;
List *act_subtree_root_List_ptr,*new_subtree_root_List_ptr,*help_ptr;
QuadTreeNodeNM *subtree_root_ptr;
build_up_root_vertex(G,T);
act_subtree_root_List.clear();
new_subtree_root_List.clear();
act_subtree_root_List.pushFront(T.get_root_ptr());
act_subtree_root_List_ptr = &act_subtree_root_List;
new_subtree_root_List_ptr = &new_subtree_root_List;
while(!act_subtree_root_List_ptr->empty())
{
while(!act_subtree_root_List_ptr->empty())
{
subtree_root_ptr = act_subtree_root_List_ptr->popFrontRet();
construct_subtree(A,T,subtree_root_ptr,*new_subtree_root_List_ptr);
}
help_ptr = act_subtree_root_List_ptr;
act_subtree_root_List_ptr = new_subtree_root_List_ptr;
new_subtree_root_List_ptr = help_ptr;
}
}
void NMM::build_up_root_vertex(const Graph&G, QuadTreeNM& T)
{
node v;
T.init_tree();
T.get_root_ptr()->set_Sm_level(0);
T.get_root_ptr()->set_Sm_downleftcorner(down_left_corner);
T.get_root_ptr()->set_Sm_boxlength(boxlength);
T.get_root_ptr()->set_particlenumber_in_subtree(G.numberOfNodes());
forall_nodes(v,G)
T.get_root_ptr()->pushBack_contained_nodes(v);
}
void NMM::construct_subtree(
NodeArray& A,
QuadTreeNM& T,
QuadTreeNodeNM *subtree_root_ptr,
List& new_subtree_root_List)
{
int n = subtree_root_ptr->get_particlenumber_in_subtree();
int subtree_depth = static_cast(max(1.0,floor(Math::log4(n))-2.0));
int maxindex=1;
for(int i=1; i<=subtree_depth; i++)
maxindex *= 2;
double subtree_min_boxlength = subtree_root_ptr->get_Sm_boxlength()/maxindex;
if(subtree_min_boxlength >= MIN_BOX_LENGTH)
{
Array2D leaf_ptr(0,maxindex-1,0,maxindex-1);
T.set_act_ptr(subtree_root_ptr);
if (find_smallest_quad(A,T)) //not all nodes have the same position
{
construct_complete_subtree(T,subtree_depth,leaf_ptr,0,0,0);
set_contained_nodes_for_leaves(A,subtree_root_ptr,leaf_ptr,maxindex);
T.set_act_ptr(subtree_root_ptr);
set_particlenumber_in_subtree_entries(T);
T.set_act_ptr(subtree_root_ptr);
construct_reduced_subtree(A,T,new_subtree_root_List);
}
}
}
void NMM::construct_complete_subtree(
QuadTreeNM& T,
int subtree_depth,
Array2D& leaf_ptr,
int act_depth,
int act_x_index,
int act_y_index)
{
if(act_depth < subtree_depth)
{
T.create_new_lt_child();
T.create_new_rt_child();
T.create_new_lb_child();
T.create_new_rb_child();
T.go_to_lt_child();
construct_complete_subtree(T,subtree_depth,leaf_ptr,act_depth+1,2*act_x_index,
2*act_y_index+1);
T.go_to_father();
T.go_to_rt_child();
construct_complete_subtree(T,subtree_depth,leaf_ptr,act_depth+1,2*act_x_index+1,
2*act_y_index+1);
T.go_to_father();
T.go_to_lb_child();
construct_complete_subtree(T,subtree_depth,leaf_ptr,act_depth+1,2*act_x_index,
2*act_y_index);
T.go_to_father();
T.go_to_rb_child();
construct_complete_subtree(T,subtree_depth,leaf_ptr,act_depth+1,2*act_x_index+1,
2*act_y_index);
T.go_to_father();
}
else if (act_depth == subtree_depth)
{
leaf_ptr(act_x_index,act_y_index) = T.get_act_ptr();
}
else
cout<<"Error NMM::construct_complete_subtree()"< &A,
QuadTreeNodeNM* subtree_root_ptr,
Array2D &leaf_ptr,
int maxindex)
{
node v;
QuadTreeNodeNM* act_ptr;
double xcoord,ycoord;
int x_index,y_index;
double minboxlength = subtree_root_ptr->get_Sm_boxlength()/maxindex;
while(!subtree_root_ptr->contained_nodes_empty())
{
v = subtree_root_ptr->pop_contained_nodes();
xcoord = A[v].get_x()-subtree_root_ptr->get_Sm_downleftcorner().m_x;
ycoord = A[v].get_y()-subtree_root_ptr->get_Sm_downleftcorner().m_y;;
x_index = int(xcoord/minboxlength);
y_index = int(ycoord/minboxlength);
act_ptr = leaf_ptr(x_index,y_index);
act_ptr->pushBack_contained_nodes(v);
act_ptr->set_particlenumber_in_subtree(act_ptr->get_particlenumber_in_subtree()+1);
}
}
void NMM::set_particlenumber_in_subtree_entries(QuadTreeNM& T)
{
int child_nr;
if(!T.get_act_ptr()->is_leaf())
{//if
T.get_act_ptr()->set_particlenumber_in_subtree(0);
if (T.get_act_ptr()->child_lt_exists())
{
T.go_to_lt_child();
set_particlenumber_in_subtree_entries(T);
T.go_to_father();
child_nr = T.get_act_ptr()->get_child_lt_ptr()->get_particlenumber_in_subtree();
T.get_act_ptr()->set_particlenumber_in_subtree(child_nr + T.get_act_ptr()->
get_particlenumber_in_subtree());
}
if (T.get_act_ptr()->child_rt_exists())
{
T.go_to_rt_child();
set_particlenumber_in_subtree_entries(T);
T.go_to_father();
child_nr = T.get_act_ptr()->get_child_rt_ptr()->get_particlenumber_in_subtree();
T.get_act_ptr()->set_particlenumber_in_subtree(child_nr + T.get_act_ptr()->
get_particlenumber_in_subtree());
}
if (T.get_act_ptr()->child_lb_exists())
{
T.go_to_lb_child();
set_particlenumber_in_subtree_entries(T);
T.go_to_father();
child_nr = T.get_act_ptr()->get_child_lb_ptr()->get_particlenumber_in_subtree();
T.get_act_ptr()->set_particlenumber_in_subtree(child_nr + T.get_act_ptr()->
get_particlenumber_in_subtree());
}
if (T.get_act_ptr()->child_rb_exists())
{
T.go_to_rb_child();
set_particlenumber_in_subtree_entries(T);
T.go_to_father();
child_nr = T.get_act_ptr()->get_child_rb_ptr()->get_particlenumber_in_subtree();
T.get_act_ptr()->set_particlenumber_in_subtree(child_nr + T.get_act_ptr()->
get_particlenumber_in_subtree());
}
}//if
}
void NMM::construct_reduced_subtree(
NodeArray& A,
QuadTreeNM& T,
List& new_subtree_root_List)
{
do
{
QuadTreeNodeNM* act_ptr = T.get_act_ptr();
delete_empty_subtrees(T);
T.set_act_ptr(act_ptr);
}
while(check_and_delete_degenerated_node(T)== true) ;
if(!T.get_act_ptr()->is_leaf() && T.get_act_ptr()->get_particlenumber_in_subtree()
<= particles_in_leaves())
{
delete_sparse_subtree(T,T.get_act_ptr());
}
//push leaves that contain many particles
if(T.get_act_ptr()->is_leaf() && T.get_act_ptr()->
get_particlenumber_in_subtree() > particles_in_leaves())
new_subtree_root_List.pushBack(T.get_act_ptr());
//find smallest quad for leaves of T
else if(T.get_act_ptr()->is_leaf() && T.get_act_ptr()->
get_particlenumber_in_subtree() <= particles_in_leaves())
find_smallest_quad(A,T);
//recursive calls
else if(!T.get_act_ptr()->is_leaf())
{//else
if(T.get_act_ptr()->child_lt_exists())
{
T.go_to_lt_child();
construct_reduced_subtree(A,T,new_subtree_root_List);
T.go_to_father();
}
if(T.get_act_ptr()->child_rt_exists())
{
T.go_to_rt_child();
construct_reduced_subtree(A,T,new_subtree_root_List);
T.go_to_father();
}
if(T.get_act_ptr()->child_lb_exists())
{
T.go_to_lb_child();
construct_reduced_subtree(A,T,new_subtree_root_List);
T.go_to_father();
}
if(T.get_act_ptr()->child_rb_exists())
{
T.go_to_rb_child();
construct_reduced_subtree(A,T,new_subtree_root_List);
T.go_to_father();
}
}//else
}
void NMM::delete_empty_subtrees(QuadTreeNM& T)
{
int child_part_nr;
QuadTreeNodeNM* act_ptr = T.get_act_ptr();
if(act_ptr->child_lt_exists())
{
child_part_nr = act_ptr->get_child_lt_ptr()->get_particlenumber_in_subtree();
if(child_part_nr == 0)
{
T.delete_tree(act_ptr->get_child_lt_ptr());
act_ptr->set_child_lt_ptr(NULL);
}
}
if(act_ptr->child_rt_exists())
{
child_part_nr = act_ptr->get_child_rt_ptr()->get_particlenumber_in_subtree();
if(child_part_nr == 0)
{
T.delete_tree(act_ptr->get_child_rt_ptr());
act_ptr->set_child_rt_ptr(NULL);
}
}
if(act_ptr->child_lb_exists())
{
child_part_nr = act_ptr->get_child_lb_ptr()->get_particlenumber_in_subtree();
if(child_part_nr == 0)
{
T.delete_tree(act_ptr->get_child_lb_ptr());
act_ptr->set_child_lb_ptr(NULL);
}
}
if(act_ptr->child_rb_exists())
{
child_part_nr = act_ptr->get_child_rb_ptr()->get_particlenumber_in_subtree();
if(child_part_nr == 0)
{
T.delete_tree(act_ptr->get_child_rb_ptr());
act_ptr->set_child_rb_ptr(NULL);
}
}
}
bool NMM::check_and_delete_degenerated_node(QuadTreeNM& T)
{
QuadTreeNodeNM* delete_ptr;
QuadTreeNodeNM* father_ptr;
QuadTreeNodeNM* child_ptr;
bool lt_child = T.get_act_ptr()->child_lt_exists();
bool rt_child = T.get_act_ptr()->child_rt_exists();
bool lb_child = T.get_act_ptr()->child_lb_exists();
bool rb_child = T.get_act_ptr()->child_rb_exists();
bool is_degenerated = false;
if(lt_child && !rt_child && !lb_child && !rb_child)
{//if1
is_degenerated = true;
delete_ptr = T.get_act_ptr();
child_ptr = T.get_act_ptr()->get_child_lt_ptr();
if(T.get_act_ptr() == T.get_root_ptr())//special case
{
T.set_root_ptr(child_ptr);
T.set_act_ptr(T.get_root_ptr());
}
else//usual case
{
father_ptr = T.get_act_ptr()->get_father_ptr();
child_ptr->set_father_ptr(father_ptr);
if(father_ptr->get_child_lt_ptr() == T.get_act_ptr())
father_ptr->set_child_lt_ptr(child_ptr);
else if(father_ptr->get_child_rt_ptr() == T.get_act_ptr())
father_ptr->set_child_rt_ptr(child_ptr);
else if(father_ptr->get_child_lb_ptr() == T.get_act_ptr())
father_ptr->set_child_lb_ptr(child_ptr);
else if(father_ptr->get_child_rb_ptr() == T.get_act_ptr())
father_ptr->set_child_rb_ptr(child_ptr);
else
cout<<"Error NMM::delete_degenerated_node"<get_child_rt_ptr();
if(T.get_act_ptr() == T.get_root_ptr())//special case
{
T.set_root_ptr(child_ptr);
T.set_act_ptr(T.get_root_ptr());
}
else//usual case
{
father_ptr = T.get_act_ptr()->get_father_ptr();
child_ptr->set_father_ptr(father_ptr);
if(father_ptr->get_child_lt_ptr() == T.get_act_ptr())
father_ptr->set_child_lt_ptr(child_ptr);
else if(father_ptr->get_child_rt_ptr() == T.get_act_ptr())
father_ptr->set_child_rt_ptr(child_ptr);
else if(father_ptr->get_child_lb_ptr() == T.get_act_ptr())
father_ptr->set_child_lb_ptr(child_ptr);
else if(father_ptr->get_child_rb_ptr() == T.get_act_ptr())
father_ptr->set_child_rb_ptr(child_ptr);
else
cout<<"Error NMM::delete_degenerated_node"<get_child_lb_ptr();
if(T.get_act_ptr() == T.get_root_ptr())//special case
{
T.set_root_ptr(child_ptr);
T.set_act_ptr(T.get_root_ptr());
}
else//usual case
{
father_ptr = T.get_act_ptr()->get_father_ptr();
child_ptr->set_father_ptr(father_ptr);
if(father_ptr->get_child_lt_ptr() == T.get_act_ptr())
father_ptr->set_child_lt_ptr(child_ptr);
else if(father_ptr->get_child_rt_ptr() == T.get_act_ptr())
father_ptr->set_child_rt_ptr(child_ptr);
else if(father_ptr->get_child_lb_ptr() == T.get_act_ptr())
father_ptr->set_child_lb_ptr(child_ptr);
else if(father_ptr->get_child_rb_ptr() == T.get_act_ptr())
father_ptr->set_child_rb_ptr(child_ptr);
else
cout<<"Error NMM::delete_degenerated_node"<get_child_rb_ptr();
if(T.get_act_ptr() == T.get_root_ptr())//special case
{
T.set_root_ptr(child_ptr);
T.set_act_ptr(T.get_root_ptr());
}
else//usual case
{
father_ptr = T.get_act_ptr()->get_father_ptr();
child_ptr->set_father_ptr(father_ptr);
if(father_ptr->get_child_lt_ptr() == T.get_act_ptr())
father_ptr->set_child_lt_ptr(child_ptr);
else if(father_ptr->get_child_rt_ptr() == T.get_act_ptr())
father_ptr->set_child_rt_ptr(child_ptr);
else if(father_ptr->get_child_lb_ptr() == T.get_act_ptr())
father_ptr->set_child_lb_ptr(child_ptr);
else if(father_ptr->get_child_rb_ptr() == T.get_act_ptr())
father_ptr->set_child_rb_ptr(child_ptr);
else
cout<<"Error NMM::delete_degenerated_node"<child_lt_exists())
{
T.delete_tree(new_leaf_ptr->get_child_lt_ptr());
new_leaf_ptr->set_child_lt_ptr(NULL);
}
if(new_leaf_ptr->child_rt_exists())
{
T.delete_tree(new_leaf_ptr->get_child_rt_ptr());
new_leaf_ptr->set_child_rt_ptr(NULL);
}
if(new_leaf_ptr->child_lb_exists())
{
T.delete_tree(new_leaf_ptr->get_child_lb_ptr());
new_leaf_ptr->set_child_lb_ptr(NULL);
}
if(new_leaf_ptr->child_rb_exists())
{
T.delete_tree(new_leaf_ptr->get_child_rb_ptr());
new_leaf_ptr->set_child_rb_ptr(NULL);
}
}
void NMM::collect_contained_nodes(QuadTreeNM& T, QuadTreeNodeNM* new_leaf_ptr)
{
if(T.get_act_ptr()->is_leaf())
while(!T.get_act_ptr()->contained_nodes_empty())
new_leaf_ptr->pushBack_contained_nodes(T.get_act_ptr()->pop_contained_nodes());
else if(T.get_act_ptr()->child_lt_exists())
{
T.go_to_lt_child();
collect_contained_nodes(T,new_leaf_ptr);
T.go_to_father();
}
if(T.get_act_ptr()->child_rt_exists())
{
T.go_to_rt_child();
collect_contained_nodes(T,new_leaf_ptr);
T.go_to_father();
}
if(T.get_act_ptr()->child_lb_exists())
{
T.go_to_lb_child();
collect_contained_nodes(T,new_leaf_ptr);
T.go_to_father();
}
if(T.get_act_ptr()->child_rb_exists())
{
T.go_to_rb_child();
collect_contained_nodes(T,new_leaf_ptr);
T.go_to_father();
}
}
bool NMM::find_smallest_quad(NodeArray& A, QuadTreeNM& T)
{
OGDF_ASSERT(!T.get_act_ptr()->contained_nodes_empty());
//if(T.get_act_ptr()->contained_nodes_empty())
// cout<<"Error NMM :: find_smallest_quad()"<L;
T.get_act_ptr()->get_contained_nodes(L);
node v = L.popFrontRet();
double x_min = A[v].get_x();
double x_max = x_min;
double y_min = A[v].get_y();
double y_max = y_min;
while(! L.empty())
{
v = L.popFrontRet();
if(A[v].get_x() < x_min)
x_min = A[v].get_x();
if(A[v].get_x() > x_max)
x_max = A[v].get_x();
if(A[v].get_y() < y_min)
y_min = A[v].get_y();
if(A[v].get_y() > y_max)
y_max = A[v].get_y();
}
if(x_min != x_max || y_min != y_max) //nodes are not all at the same position
{
if(find_sm_cell() == FMMMLayout::scfIteratively)
find_small_cell_iteratively(T.get_act_ptr(),x_min,x_max,y_min,y_max);
else //find_sm_cell == FMMMLayout::scfAluru
find_small_cell_iteratively(T.get_act_ptr(),x_min,x_max,y_min,y_max);
return true;
}
else
return false;
//}//else
}
// ********Functions needed for subtree by subtree tree construction(END)************
void NMM::find_small_cell_iteratively(
QuadTreeNodeNM* act_ptr,
double x_min,
double x_max,
double y_min,
double y_max)
{
int new_level;
double new_boxlength;
DPoint new_dlc;
bool Sm_cell_found = false;
while ( !Sm_cell_found && ((x_max-x_min >=MIN_BOX_LENGTH) ||
(y_max-y_min >=MIN_BOX_LENGTH)) )
{
if(in_lt_quad(act_ptr,x_min,x_max,y_min,y_max))
{
new_level = act_ptr->get_Sm_level()+1;
new_boxlength = act_ptr->get_Sm_boxlength()/2;
new_dlc.m_x = act_ptr->get_Sm_downleftcorner().m_x;
new_dlc.m_y = act_ptr->get_Sm_downleftcorner().m_y+new_boxlength;
act_ptr->set_Sm_level(new_level);
act_ptr->set_Sm_boxlength(new_boxlength);
act_ptr->set_Sm_downleftcorner(new_dlc);
}
else if(in_rt_quad(act_ptr,x_min,x_max,y_min,y_max))
{
new_level = act_ptr->get_Sm_level()+1;
new_boxlength = act_ptr->get_Sm_boxlength()/2;
new_dlc.m_x = act_ptr->get_Sm_downleftcorner().m_x+new_boxlength;
new_dlc.m_y = act_ptr->get_Sm_downleftcorner().m_y+new_boxlength;
act_ptr->set_Sm_level(new_level);
act_ptr->set_Sm_boxlength(new_boxlength);
act_ptr->set_Sm_downleftcorner(new_dlc);
}
else if(in_lb_quad(act_ptr,x_min,x_max,y_min,y_max))
{
new_level = act_ptr->get_Sm_level()+1;
new_boxlength = act_ptr->get_Sm_boxlength()/2;
act_ptr->set_Sm_level(new_level);
act_ptr->set_Sm_boxlength(new_boxlength);
}
else if(in_rb_quad(act_ptr,x_min,x_max,y_min,y_max))
{
new_level = act_ptr->get_Sm_level()+1;
new_boxlength = act_ptr->get_Sm_boxlength()/2;
new_dlc.m_x = act_ptr->get_Sm_downleftcorner().m_x+new_boxlength;
new_dlc.m_y = act_ptr->get_Sm_downleftcorner().m_y;
act_ptr->set_Sm_level(new_level);
act_ptr->set_Sm_boxlength(new_boxlength);
act_ptr->set_Sm_downleftcorner(new_dlc);
}
else Sm_cell_found = true;
}
}
void NMM::find_small_cell_by_formula(
QuadTreeNodeNM* act_ptr,
double x_min,
double x_max,
double y_min,
double y_max)
{
numexcept N;
int level_offset = act_ptr->get_Sm_level();
max_power_of_2_index = 30;//up to this level standard integer arithmetic is used
// DPoint nullpoint (0,0);
IPoint Sm_position;
double Sm_dlc_x_coord,Sm_dlc_y_coord;
double Sm_boxlength;
int Sm_level;
DPoint Sm_downleftcorner;
int j_x = max_power_of_2_index+1;
int j_y = max_power_of_2_index+1;
bool rectangle_is_horizontal_line = false;
bool rectangle_is_vertical_line = false;
bool rectangle_is_point = false;
//shift boundaries to the origin for easy calculations
double x_min_old = x_min;
double x_max_old = x_max;
double y_min_old = y_min;
double y_max_old = y_max;
Sm_boxlength = act_ptr->get_Sm_boxlength();
Sm_dlc_x_coord = act_ptr->get_Sm_downleftcorner().m_x;
Sm_dlc_y_coord = act_ptr->get_Sm_downleftcorner().m_y;
x_min -= Sm_dlc_x_coord;
x_max -= Sm_dlc_x_coord;
y_min -= Sm_dlc_y_coord;
y_max -= Sm_dlc_y_coord;
//check if iterative way has to be used
if (x_min == x_max && y_min == y_max)
rectangle_is_point = true;
else if(x_min == x_max && y_min != y_max)
rectangle_is_vertical_line = true;
else //x_min != x_max
j_x = static_cast(ceil(Math::log2(Sm_boxlength/(x_max-x_min))));
if(x_min != x_max && y_min == y_max)
rectangle_is_horizontal_line = true;
else //y_min != y_max
j_y = static_cast(ceil(Math::log2(Sm_boxlength/(y_max-y_min))));
if(rectangle_is_point)
{
;//keep the old values
}
else if ( !N.nearly_equal((x_min_old - x_max_old),(x_min-x_max)) ||
!N.nearly_equal((y_min_old - y_max_old),(y_min-y_max)) ||
x_min/Sm_boxlength < MIN_BOX_LENGTH || x_max/Sm_boxlength < MIN_BOX_LENGTH ||
y_min/Sm_boxlength < MIN_BOX_LENGTH || y_max/Sm_boxlength < MIN_BOX_LENGTH )
find_small_cell_iteratively(act_ptr,x_min_old,x_max_old,y_min_old,y_max_old);
else if ( ((j_x > max_power_of_2_index) && (j_y > max_power_of_2_index)) ||
((j_x > max_power_of_2_index) && !rectangle_is_vertical_line) ||
((j_y > max_power_of_2_index) && !rectangle_is_horizontal_line) )
find_small_cell_iteratively(act_ptr,x_min_old,x_max_old,y_min_old,y_max_old);
else //idea of Aluru et al.
{//else
int k,a1,a2,A,j_minus_k;
double h1,h2;
int Sm_x_level = 0, Sm_y_level = 0;
int Sm_x_position = 0, Sm_y_position = 0;
if(x_min != x_max)
{//if1
//calculate Sm_x_level and Sm_x_position
a1 = static_cast(ceil((x_min/Sm_boxlength)*power_of_two(j_x)));
a2 = static_cast(floor((x_max/Sm_boxlength)*power_of_two(j_x)));
h1 = (Sm_boxlength/power_of_two(j_x))* a1;
h2 = (Sm_boxlength/power_of_two(j_x))* a2;
//special cases: two tangents or left tangent and righ cutline
if(((h1 == x_min)&&(h2 == x_max)) || ((h1 == x_min) && (h2 != x_max)) )
A = a2;
else if (a1 == a2) //only one cutline
A = a1;
else //two cutlines or a right tangent and a left cutline (usual case)
{
if((a1 % 2) == 0)
A = a1;
else
A = a2;
}
j_minus_k = static_cast(Math::log2(1+(A ^ (A-1)))-1);
k = j_x - j_minus_k;
Sm_x_level = k-1;
Sm_x_position = a1/ power_of_two(j_x - Sm_x_level);
}//if1
if(y_min != y_max)
{//if2
//calculate Sm_y_level and Sm_y_position
a1 = static_cast(ceil((y_min/Sm_boxlength)*power_of_two(j_y)));
a2 = static_cast(floor((y_max/Sm_boxlength)*power_of_two(j_y)));
h1 = (Sm_boxlength/power_of_two(j_y))* a1;
h2 = (Sm_boxlength/power_of_two(j_y))* a2;
//special cases: two tangents or bottom tangent and top cutline
if(((h1 == y_min)&&(h2 == y_max)) || ((h1 == y_min) && (h2 != y_max)) )
A = a2;
else if (a1 == a2) //only one cutline
A = a1;
else //two cutlines or a top tangent and a bottom cutline (usual case)
{
if((a1 % 2) == 0)
A = a1;
else
A = a2;
}
j_minus_k = static_cast(Math::log2(1+(A ^ (A-1)))-1);
k = j_y - j_minus_k;
Sm_y_level = k-1;
Sm_y_position = a1/ power_of_two(j_y - Sm_y_level);
}//if2
if((x_min != x_max) &&(y_min != y_max))//a box with area > 0
{//if3
if (Sm_x_level == Sm_y_level)
{
Sm_level = Sm_x_level;
Sm_position.m_x = Sm_x_position;
Sm_position.m_y = Sm_y_position;
}
else if (Sm_x_level < Sm_y_level)
{
Sm_level = Sm_x_level;
Sm_position.m_x = Sm_x_position;
Sm_position.m_y = Sm_y_position/power_of_two(Sm_y_level-Sm_x_level);
}
else //Sm_x_level > Sm_y_level
{
Sm_level = Sm_y_level;
Sm_position.m_x = Sm_x_position/power_of_two(Sm_x_level-Sm_y_level);
Sm_position.m_y = Sm_y_position;
}
}//if3
else if(x_min == x_max) //a vertical line
{//if4
Sm_level = Sm_y_level;
Sm_position.m_x = static_cast (floor((x_min*power_of_two(Sm_level))/
Sm_boxlength));
Sm_position.m_y = Sm_y_position;
}//if4
else //y_min == y_max (a horizontal line)
{//if5
Sm_level = Sm_x_level;
Sm_position.m_x = Sm_x_position;
Sm_position.m_y = static_cast (floor((y_min*power_of_two(Sm_level))/
Sm_boxlength));
}//if5
Sm_boxlength = Sm_boxlength/power_of_two(Sm_level);
Sm_downleftcorner.m_x = Sm_dlc_x_coord + Sm_boxlength * Sm_position.m_x;
Sm_downleftcorner.m_y = Sm_dlc_y_coord + Sm_boxlength * Sm_position.m_y;
act_ptr->set_Sm_level(Sm_level+level_offset);
act_ptr->set_Sm_boxlength(Sm_boxlength);
act_ptr->set_Sm_downleftcorner(Sm_downleftcorner);
}//else
}
inline void NMM::delete_red_quad_tree_and_count_treenodes(QuadTreeNM& T)
{
T.delete_tree(T.get_root_ptr());
}
inline void NMM::form_multipole_expansions(
NodeArray& A,
QuadTreeNM& T,
List& quad_tree_leaves)
{
T.set_act_ptr(T.get_root_ptr());
form_multipole_expansion_of_subtree(A,T,quad_tree_leaves);
}
void NMM::form_multipole_expansion_of_subtree(
NodeArray& A,
QuadTreeNM& T,
List& quad_tree_leaves)
{
init_expansion_Lists(T.get_act_ptr());
set_center(T.get_act_ptr());
if(T.get_act_ptr()->is_leaf()) //form expansions for leaf nodes
{//if
quad_tree_leaves.pushBack(T.get_act_ptr());
form_multipole_expansion_of_leaf_node(A,T.get_act_ptr());
}//if
else //rekursive calls and add shifted expansions
{//else
if(T.get_act_ptr()->child_lt_exists())
{
T.go_to_lt_child();
form_multipole_expansion_of_subtree(A,T,quad_tree_leaves);
add_shifted_expansion_to_father_expansion(T.get_act_ptr());
T.go_to_father();
}
if(T.get_act_ptr()->child_rt_exists())
{
T.go_to_rt_child();
form_multipole_expansion_of_subtree(A,T,quad_tree_leaves);
add_shifted_expansion_to_father_expansion(T.get_act_ptr());
T.go_to_father();
}
if(T.get_act_ptr()->child_lb_exists())
{
T.go_to_lb_child();
form_multipole_expansion_of_subtree(A,T,quad_tree_leaves);
add_shifted_expansion_to_father_expansion(T.get_act_ptr());
T.go_to_father();
}
if(T.get_act_ptr()->child_rb_exists())
{
T.go_to_rb_child();
form_multipole_expansion_of_subtree(A,T,quad_tree_leaves);
add_shifted_expansion_to_father_expansion(T.get_act_ptr());
T.go_to_father();
}
}//else
}
inline void NMM::init_expansion_Lists(QuadTreeNodeNM* act_ptr)
{
int i;
Array > nulList (precision()+1);
for (i = 0;i<=precision();i++)
nulList[i] = 0;
act_ptr->set_multipole_exp(nulList,precision());
act_ptr->set_locale_exp(nulList,precision());
}
void NMM::set_center(QuadTreeNodeNM* act_ptr)
{
const int BILLION = 1000000000;
DPoint Sm_downleftcorner = act_ptr->get_Sm_downleftcorner();
double Sm_boxlength = act_ptr->get_Sm_boxlength();
double boxcenter_x_coord,boxcenter_y_coord;
DPoint Sm_dlc;
double rand_y;
boxcenter_x_coord = Sm_downleftcorner.m_x + Sm_boxlength * 0.5;
boxcenter_y_coord = Sm_downleftcorner.m_y + Sm_boxlength * 0.5;
//for use of complex logarithm: waggle the y-coordinates a little bit
//such that the new center is really inside the actual box and near the exact center
rand_y = double(randomNumber(1,BILLION)+1)/(BILLION+2);//rand number in (0,1)
boxcenter_y_coord = boxcenter_y_coord + 0.001 * Sm_boxlength * rand_y;
complex boxcenter(boxcenter_x_coord,boxcenter_y_coord);
act_ptr->set_Sm_center(boxcenter);
}
void NMM::form_multipole_expansion_of_leaf_node(
NodeArray& A,
QuadTreeNodeNM* act_ptr)
{
int k;
complex Q (0,0);
complex z_0 = act_ptr->get_Sm_center();//center of actual box
complex nullpoint (0,0);
Array > coef (precision()+1);
complex z_v_minus_z_0_over_k;
List nodes_in_box;
int i;
ListIterator v_it;
act_ptr->get_contained_nodes(nodes_in_box);
for(v_it = nodes_in_box.begin();v_it.valid();++v_it)
Q += 1;
coef[0] = Q;
for(i = 1; i<=precision();i++)
coef[i] = nullpoint;
for(v_it = nodes_in_box.begin();v_it.valid();++v_it)
{
complex z_v (A[*v_it].get_x(),A[*v_it].get_y());
z_v_minus_z_0_over_k = z_v - z_0;
for(k=1;k<=precision();k++)
{
coef[k] += ((double(-1))*z_v_minus_z_0_over_k)/double(k);
z_v_minus_z_0_over_k *= z_v - z_0;
}
}
act_ptr->replace_multipole_exp(coef,precision());
}
void NMM::add_shifted_expansion_to_father_expansion(QuadTreeNodeNM* act_ptr)
{
QuadTreeNodeNM* father_ptr = act_ptr->get_father_ptr();
complex sum;
complex z_0,z_1;
Array > z_0_minus_z_1_over (precision()+1);
z_1 = father_ptr->get_Sm_center();
z_0 = act_ptr->get_Sm_center();
father_ptr->get_multipole_exp()[0] += act_ptr->get_multipole_exp()[0];
//init z_0_minus_z_1_over
z_0_minus_z_1_over[0] = 1;
for(int i = 1; i<= precision(); i++)
z_0_minus_z_1_over[i] = z_0_minus_z_1_over[i-1] * (z_0 - z_1);
for(int k=1; k<=precision(); k++)
{
sum = (act_ptr->get_multipole_exp()[0]*(double(-1))*z_0_minus_z_1_over[k])/
double(k) ;
for(int s=1; s<=k; s++)
sum += act_ptr->get_multipole_exp()[s]*z_0_minus_z_1_over[k-s]* binko(k-1,s-1);
father_ptr->get_multipole_exp()[k] += sum;
}
}
void NMM::calculate_local_expansions_and_WSPRLS(
NodeArray&A,
QuadTreeNodeNM* act_node_ptr)
{
List I,L,L2,E,D1,D2,M;
QuadTreeNodeNM *father_ptr = 0, *selected_node_ptr;
ListIterator ptr_it;
//Step 0: Initializations
if(! act_node_ptr->is_root())
father_ptr = act_node_ptr->get_father_ptr();
I.clear();L.clear();L2.clear();D1.clear();D2.clear();M.clear();
//Step 1: calculate Lists I (min. ill sep. set), L (interaction List of well sep.
//nodes , they are used to form the Local Expansions from the multipole expansions),
//L2 (non bordering leaves that have a larger or equal Sm-cell and are ill separated;
//empty if the actual node is a leaf)
//calculate List D1(bordering leaves that have a larger or equal Sm-cell and are
//ill separated) and D2 (non bordering leaves that have a larger or equal Sm-cell and
//are ill separated;empty if the actual node is an interior node)
//special case: act_node is the root of T
if (act_node_ptr->is_root())
{//if
E.clear();
if(act_node_ptr->child_lt_exists())
E.pushBack(act_node_ptr->get_child_lt_ptr());
if(act_node_ptr->child_rt_exists())
E.pushBack(act_node_ptr->get_child_rt_ptr());
if(act_node_ptr->child_lb_exists())
E.pushBack(act_node_ptr->get_child_lb_ptr());
if(act_node_ptr->child_rb_exists())
E.pushBack(act_node_ptr->get_child_rb_ptr());
}//if
//usual case: act_node is an interior node of T
else
{
father_ptr->get_D1(E); //bordering leaves of father
father_ptr->get_I(I); //min ill sep. nodes of father
for(ptr_it = I.begin();ptr_it.valid();++ptr_it)
E.pushBack(*ptr_it);
I.clear();
}
while (!E.empty())
{//while
selected_node_ptr = E.popFrontRet();
if (well_separated(act_node_ptr,selected_node_ptr))
L.pushBack(selected_node_ptr);
else if (act_node_ptr->get_Sm_level() < selected_node_ptr->get_Sm_level())
I.pushBack(selected_node_ptr);
else if(!selected_node_ptr->is_leaf())
{
if(selected_node_ptr->child_lt_exists())
E.pushBack(selected_node_ptr->get_child_lt_ptr());
if(selected_node_ptr->child_rt_exists())
E.pushBack(selected_node_ptr->get_child_rt_ptr());
if(selected_node_ptr->child_lb_exists())
E.pushBack(selected_node_ptr->get_child_lb_ptr());
if(selected_node_ptr->child_rb_exists())
E.pushBack(selected_node_ptr->get_child_rb_ptr());
}
else if(bordering(act_node_ptr,selected_node_ptr))
D1.pushBack(selected_node_ptr);
else if( (selected_node_ptr != act_node_ptr)&&(act_node_ptr->is_leaf()))
D2.pushBack(selected_node_ptr); //direct calculation (no errors produced)
else if((selected_node_ptr != act_node_ptr)&&!(act_node_ptr->is_leaf()))
L2.pushBack(selected_node_ptr);
}//while
act_node_ptr->set_I(I);
act_node_ptr->set_D1(D1);
act_node_ptr->set_D2(D2);
//Step 2: add local expansions from father(act_node_ptr) and calculate locale
//expansions for all nodes in L
if(!act_node_ptr->is_root())
add_shifted_local_exp_of_parent(act_node_ptr);
for(ptr_it = L.begin();ptr_it.valid();++ptr_it)
add_local_expansion(*ptr_it,act_node_ptr);
//Step 3: calculate locale expansions for all nodes in D2 (simpler than in Step 2)
for(ptr_it = L2.begin();ptr_it.valid();++ptr_it)
add_local_expansion_of_leaf(A,*ptr_it,act_node_ptr);
//Step 4: recursive calls if act_node is not a leaf
if(!act_node_ptr->is_leaf())
{
if(act_node_ptr->child_lt_exists())
calculate_local_expansions_and_WSPRLS(A,act_node_ptr->get_child_lt_ptr());
if(act_node_ptr->child_rt_exists())
calculate_local_expansions_and_WSPRLS(A,act_node_ptr->get_child_rt_ptr());
if(act_node_ptr->child_lb_exists())
calculate_local_expansions_and_WSPRLS(A,act_node_ptr->get_child_lb_ptr());
if(act_node_ptr->child_rb_exists())
calculate_local_expansions_and_WSPRLS(A,act_node_ptr->get_child_rb_ptr());
}
//Step 5: WSPRLS(Well Separateness Preserving Refinement of leaf surroundings)
//if act_node is a leaf than calculate the list D1,D2 and M from I and D1
else // *act_node_ptr is a leaf
{//else
act_node_ptr->get_D1(D1);
act_node_ptr->get_D2(D2);
while(!I.empty())
{//while
selected_node_ptr = I.popFrontRet();
if(selected_node_ptr->is_leaf())
{
//here D1 contains larger AND smaller bordering leaves!
if(bordering(act_node_ptr,selected_node_ptr))
D1.pushBack(selected_node_ptr);
else
D2.pushBack(selected_node_ptr);
}
else //!selected_node_ptr->is_leaf()
{
if(bordering(act_node_ptr,selected_node_ptr))
{
if(selected_node_ptr->child_lt_exists())
I.pushBack(selected_node_ptr->get_child_lt_ptr());
if(selected_node_ptr->child_rt_exists())
I.pushBack(selected_node_ptr->get_child_rt_ptr());
if(selected_node_ptr->child_lb_exists())
I.pushBack(selected_node_ptr->get_child_lb_ptr());
if(selected_node_ptr->child_rb_exists())
I.pushBack(selected_node_ptr->get_child_rb_ptr());
}
else
M.pushBack(selected_node_ptr);
}
}//while
act_node_ptr->set_D1(D1);
act_node_ptr->set_D2(D2);
act_node_ptr->set_M(M);
}//else
}
bool NMM::well_separated(QuadTreeNodeNM* node_1_ptr, QuadTreeNodeNM* node_2_ptr)
{
numexcept N;
double boxlength_1 = node_1_ptr->get_Sm_boxlength();
double boxlength_2 = node_2_ptr->get_Sm_boxlength();
double x1_min,x1_max,y1_min,y1_max,x2_min,x2_max,y2_min,y2_max;
bool x_overlap,y_overlap;
if(boxlength_1 <= boxlength_2)
{
x1_min = node_1_ptr->get_Sm_downleftcorner().m_x;
x1_max = node_1_ptr->get_Sm_downleftcorner().m_x+boxlength_1;
y1_min = node_1_ptr->get_Sm_downleftcorner().m_y;
y1_max = node_1_ptr->get_Sm_downleftcorner().m_y+boxlength_1;
//blow the box up
x2_min = node_2_ptr->get_Sm_downleftcorner().m_x-boxlength_2;
x2_max = node_2_ptr->get_Sm_downleftcorner().m_x+2*boxlength_2;
y2_min = node_2_ptr->get_Sm_downleftcorner().m_y-boxlength_2;
y2_max = node_2_ptr->get_Sm_downleftcorner().m_y+2*boxlength_2;
}
else //boxlength_1 > boxlength_2
{
//blow the box up
x1_min = node_1_ptr->get_Sm_downleftcorner().m_x-boxlength_1;
x1_max = node_1_ptr->get_Sm_downleftcorner().m_x+2*boxlength_1;
y1_min = node_1_ptr->get_Sm_downleftcorner().m_y-boxlength_1;
y1_max = node_1_ptr->get_Sm_downleftcorner().m_y+2*boxlength_1;
x2_min = node_2_ptr->get_Sm_downleftcorner().m_x;
x2_max = node_2_ptr->get_Sm_downleftcorner().m_x+boxlength_2;
y2_min = node_2_ptr->get_Sm_downleftcorner().m_y;
y2_max = node_2_ptr->get_Sm_downleftcorner().m_y+boxlength_2;
}
//test if boxes overlap
if((x1_max <= x2_min)|| N.nearly_equal(x1_max,x2_min)||
(x2_max <= x1_min)|| N.nearly_equal(x2_max,x1_min))
x_overlap = false;
else
x_overlap = true;
if((y1_max <= y2_min)|| N.nearly_equal(y1_max,y2_min)||
(y2_max <= y1_min)|| N.nearly_equal(y2_max,y1_min))
y_overlap = false;
else
y_overlap = true;
if (x_overlap && y_overlap)
return false;
else
return true;
}
bool NMM::bordering(QuadTreeNodeNM* node_1_ptr,QuadTreeNodeNM* node_2_ptr)
{
numexcept N;
double boxlength_1 = node_1_ptr->get_Sm_boxlength();
double boxlength_2 = node_2_ptr->get_Sm_boxlength();
double x1_min = node_1_ptr->get_Sm_downleftcorner().m_x;
double x1_max = node_1_ptr->get_Sm_downleftcorner().m_x+boxlength_1;
double y1_min = node_1_ptr->get_Sm_downleftcorner().m_y;
double y1_max = node_1_ptr->get_Sm_downleftcorner().m_y+boxlength_1;
double x2_min = node_2_ptr->get_Sm_downleftcorner().m_x;
double x2_max = node_2_ptr->get_Sm_downleftcorner().m_x+boxlength_2;
double y2_min = node_2_ptr->get_Sm_downleftcorner().m_y;
double y2_max = node_2_ptr->get_Sm_downleftcorner().m_y+boxlength_2;
if( ( (x2_min <= x1_min || N.nearly_equal(x2_min,x1_min)) &&
(x1_max <= x2_max || N.nearly_equal(x1_max,x2_max)) &&
(y2_min <= y1_min || N.nearly_equal(y2_min,y1_min)) &&
(y1_max <= y2_max || N.nearly_equal(y1_max,y2_max)) ) ||
( (x1_min <= x2_min || N.nearly_equal(x1_min,x2_min)) &&
(x2_max <= x1_max || N.nearly_equal(x2_max,x1_max)) &&
(y1_min <= y2_min || N.nearly_equal(y1_min,y2_min)) &&
(y2_max <= y1_max || N.nearly_equal(y2_max,y1_max)) ) )
return false; //one box contains the other box(inclusive neighbours)
else
{//else
if (boxlength_1 <= boxlength_2)
{ //shift box1
if (x1_min < x2_min)
{ x1_min +=boxlength_1;x1_max +=boxlength_1; }
else if (x1_max > x2_max)
{ x1_min -=boxlength_1;x1_max -=boxlength_1; }
if (y1_min < y2_min)
{ y1_min +=boxlength_1;y1_max +=boxlength_1; }
else if (y1_max > y2_max)
{ y1_min -=boxlength_1;y1_max -=boxlength_1; }
}
else //boxlength_1 > boxlength_2
{//shift box2
if (x2_min < x1_min)
{ x2_min +=boxlength_2;x2_max +=boxlength_2; }
else if (x2_max > x1_max)
{ x2_min -=boxlength_2;x2_max -=boxlength_2; }
if (y2_min < y1_min)
{ y2_min +=boxlength_2;y2_max +=boxlength_2; }
else if (y2_max > y1_max)
{ y2_min -=boxlength_2;y2_max -=boxlength_2; }
}
if( ( (x2_min <= x1_min || N.nearly_equal(x2_min,x1_min)) &&
(x1_max <= x2_max || N.nearly_equal(x1_max,x2_max)) &&
(y2_min <= y1_min || N.nearly_equal(y2_min,y1_min)) &&
(y1_max <= y2_max || N.nearly_equal(y1_max,y2_max)) ) ||
( (x1_min <= x2_min || N.nearly_equal(x1_min,x2_min)) &&
(x2_max <= x1_max || N.nearly_equal(x2_max,x1_max)) &&
(y1_min <= y2_min || N.nearly_equal(y1_min,y2_min)) &&
(y2_max <= y1_max || N.nearly_equal(y2_max,y1_max)) ) )
return true;
else
return false;
}//else
}
void NMM::add_shifted_local_exp_of_parent(QuadTreeNodeNM* node_ptr)
{
QuadTreeNodeNM* father_ptr = node_ptr->get_father_ptr();
complex z_0 = father_ptr->get_Sm_center();
complex z_1 = node_ptr->get_Sm_center();
Array > z_1_minus_z_0_over (precision()+1);
//init z_1_minus_z_0_over
z_1_minus_z_0_over[0] = 1;
for(int i = 1; i<= precision(); i++)
z_1_minus_z_0_over[i] = z_1_minus_z_0_over[i-1] * (z_1 - z_0);
for(int l = 0; l <= precision();l++)
{
complex sum (0,0);
for(int k = l;k<=precision();k++)
sum += binko(k,l)*father_ptr->get_local_exp()[k]*z_1_minus_z_0_over[k-l];
node_ptr->get_local_exp()[l] += sum;
}
}
void NMM::add_local_expansion(QuadTreeNodeNM* ptr_0, QuadTreeNodeNM* ptr_1)
{
complex z_0 = ptr_0->get_Sm_center();
complex z_1 = ptr_1->get_Sm_center();
complex sum, z_error;
complex factor;
complex z_1_minus_z_0_over_k;
complex z_1_minus_z_0_over_s;
complex pow_minus_1_s_plus_1;
complex pow_minus_1_s;
//Error-Handling for complex logarithm
if ((std::real(z_1-z_0) <=0) && (std::imag(z_1-z_0) == 0)) //no cont. compl. log fct exists !!!
{
z_error = log(z_1 -z_0 + 0.0000001);
sum = ptr_0->get_multipole_exp()[0] * z_error;
}
else
sum = ptr_0->get_multipole_exp()[0]* log(z_1-z_0);
z_1_minus_z_0_over_k = z_1 - z_0;
for(int k = 1; k<=precision(); k++)
{
sum += ptr_0->get_multipole_exp()[k]/z_1_minus_z_0_over_k;
z_1_minus_z_0_over_k *= z_1-z_0;
}
ptr_1->get_local_exp()[0] += sum;
z_1_minus_z_0_over_s = z_1 - z_0;
for (int s = 1; s <= precision(); s++)
{
pow_minus_1_s_plus_1 = (((s+1)% 2 == 0) ? 1 : -1);
pow_minus_1_s = ((pow_minus_1_s_plus_1 == double(1))? -1 : 1);
sum = pow_minus_1_s_plus_1*ptr_0->get_multipole_exp()[0]/(z_1_minus_z_0_over_s *
double(s));
factor = pow_minus_1_s/z_1_minus_z_0_over_s;
z_1_minus_z_0_over_s *= z_1-z_0;
complex sum_2 (0,0);
z_1_minus_z_0_over_k = z_1 - z_0;
for(int k=1; k<=precision(); k++)
{
sum_2 += binko(s+k-1,k-1)*ptr_0->get_multipole_exp()[k]/z_1_minus_z_0_over_k;
z_1_minus_z_0_over_k *= z_1-z_0;
}
ptr_1->get_local_exp()[s] += sum + factor* sum_2;
}
}
void NMM::add_local_expansion_of_leaf(
NodeArray&A,
QuadTreeNodeNM* ptr_0,
QuadTreeNodeNM* ptr_1)
{
List contained_nodes;
double multipole_0_of_v = 1;//only the first coefficient is not zero
complex z_1 = ptr_1->get_Sm_center();
complex z_error;
complex z_1_minus_z_0_over_s;
complex pow_minus_1_s_plus_1;
ptr_0->get_contained_nodes(contained_nodes);
forall_listiterators(node, v_it, contained_nodes)
{//forall
//set position of v as center ( (1,0,....,0) are the multipole coefficients at v)
complex z_0 (A[*v_it].get_x(),A[*v_it].get_y());
//now transform multipole_0_of_v to the locale expansion around z_1
//Error-Handling for complex logarithm
if ((std::real(z_1-z_0) <=0) && (std::imag(z_1-z_0) == 0)) //no cont. compl. log fct exists!
{
z_error = log(z_1 -z_0 + 0.0000001);
ptr_1->get_local_exp()[0] += multipole_0_of_v * z_error;
}
else
ptr_1->get_local_exp()[0] += multipole_0_of_v * log(z_1-z_0);
z_1_minus_z_0_over_s = z_1 - z_0;
for (int s = 1;s <= precision();s++)
{
pow_minus_1_s_plus_1 = (((s+1)% 2 == 0) ? 1 : -1);
ptr_1->get_local_exp()[s] += pow_minus_1_s_plus_1*multipole_0_of_v/
(z_1_minus_z_0_over_s * double(s));
z_1_minus_z_0_over_s *= z_1-z_0;
}
}//forall
}
void NMM::transform_local_exp_to_forces(
NodeArray &A,
List& quad_tree_leaves,
NodeArray& F_local_exp)
{
List contained_nodes;
complex sum;
complex complex_null (0,0);
complex z_0;
complex z_v_minus_z_0_over_k_minus_1;
DPoint force_vector;
//calculate derivative of the potential polynom (= local expansion at leaf nodes)
//and evaluate it for each node in contained_nodes()
//and transform the complex number back to the real-world, to obtain the force
forall_listiterators( QuadTreeNodeNM*, leaf_ptr_ptr,quad_tree_leaves)
{
(*leaf_ptr_ptr)->get_contained_nodes(contained_nodes);
z_0 = (*leaf_ptr_ptr)->get_Sm_center();
forall_listiterators(node, v_ptr,contained_nodes)
{
complex z_v (A[*v_ptr].get_x(),A[*v_ptr].get_y());
sum = complex_null;
z_v_minus_z_0_over_k_minus_1 = 1;
for(int k=1; k<=precision(); k++)
{
sum += double(k) * (*leaf_ptr_ptr)->get_local_exp()[k] *
z_v_minus_z_0_over_k_minus_1;
z_v_minus_z_0_over_k_minus_1 *= z_v - z_0;
}
force_vector.m_x = sum.real();
force_vector.m_y = (-1) * sum.imag();
F_local_exp[*v_ptr] = force_vector;
}
}
}
void NMM::transform_multipole_exp_to_forces(
NodeArray& A,
List& quad_tree_leaves,
NodeArray& F_multipole_exp)
{
List M;
List act_contained_nodes;
ListIterator v_ptr;
complex sum;
complex z_0;
complex z_v_minus_z_0_over_minus_k_minus_1;
DPoint force_vector;
//for each leaf u in the M-List of an actual leaf v do:
//calculate derivative of the multipole expansion function at u
//and evaluate it for each node in v.get_contained_nodes()
//and transform the complex number back to the real-world, to obtain the force
forall_listiterators(QuadTreeNodeNM*, act_leaf_ptr_ptr,quad_tree_leaves)
{
(*act_leaf_ptr_ptr)->get_contained_nodes(act_contained_nodes);
(*act_leaf_ptr_ptr)->get_M(M);
forall_listiterators(QuadTreeNodeNM*, M_node_ptr_ptr,M)
{
z_0 = (*M_node_ptr_ptr)->get_Sm_center();
forall_listiterators(node, v_ptr,act_contained_nodes)
{
complex z_v (A[*v_ptr].get_x(),A[*v_ptr].get_y());
z_v_minus_z_0_over_minus_k_minus_1 = 1.0/(z_v-z_0);
sum = (*M_node_ptr_ptr)->get_multipole_exp()[0]*
z_v_minus_z_0_over_minus_k_minus_1;
for(int k=1; k<=precision(); k++)
{
z_v_minus_z_0_over_minus_k_minus_1 /= z_v - z_0;
sum -= double(k) * (*M_node_ptr_ptr)->get_multipole_exp()[k] *
z_v_minus_z_0_over_minus_k_minus_1;
}
force_vector.m_x = sum.real();
force_vector.m_y = (-1) * sum.imag();
F_multipole_exp[*v_ptr] = F_multipole_exp[*v_ptr] + force_vector;
}
}
}
}
void NMM::calculate_neighbourcell_forces(
NodeArray& A,
List & quad_tree_leaves,
NodeArray& F_direct)
{
numexcept N;
List act_contained_nodes,neighbour_contained_nodes,non_neighbour_contained_nodes;
List neighboured_leaves;
List non_neighboured_leaves;
double act_leaf_boxlength,neighbour_leaf_boxlength;
DPoint act_leaf_dlc,neighbour_leaf_dlc;
DPoint f_rep_u_on_v;
DPoint vector_v_minus_u;
// DPoint nullpoint(0,0);
DPoint pos_u,pos_v;
double norm_v_minus_u,scalar;
int length;
node u,v;
forall_listiterators(QuadTreeNodeNM*, act_leaf_ptr,quad_tree_leaves)
{//forall
(*act_leaf_ptr)->get_contained_nodes(act_contained_nodes);
if(act_contained_nodes.size() <= particles_in_leaves())
{//if (usual case)
//Step1:calculate forces inside act_contained_nodes
length = act_contained_nodes.size();
Array numbered_nodes (length+1);
int k = 1;
forall_listiterators(node, v_ptr,act_contained_nodes)
{
numbered_nodes[k]= *v_ptr;
k++;
}
for(k = 1; kget_D1()
(*act_leaf_ptr)->get_D1(neighboured_leaves);
act_leaf_boxlength = (*act_leaf_ptr)->get_Sm_boxlength();
act_leaf_dlc = (*act_leaf_ptr)->get_Sm_downleftcorner();
forall_listiterators(QuadTreeNodeNM*, neighbour_leaf_ptr,neighboured_leaves)
{//forall2
//forget boxes that have already been looked at
neighbour_leaf_boxlength = (*neighbour_leaf_ptr)->get_Sm_boxlength();
neighbour_leaf_dlc = (*neighbour_leaf_ptr)->get_Sm_downleftcorner();
if( (act_leaf_boxlength > neighbour_leaf_boxlength) ||
(act_leaf_boxlength == neighbour_leaf_boxlength &&
act_leaf_dlc.m_x < neighbour_leaf_dlc.m_x)
|| (act_leaf_boxlength == neighbour_leaf_boxlength &&
act_leaf_dlc.m_x == neighbour_leaf_dlc.m_x &&
act_leaf_dlc.m_y < neighbour_leaf_dlc.m_y) )
{//if
(*neighbour_leaf_ptr)->get_contained_nodes(neighbour_contained_nodes);
forall_listiterators(node, v_ptr,act_contained_nodes)
forall_listiterators(node, u_ptr, neighbour_contained_nodes)
{//for
pos_u = A[*u_ptr].get_position();
pos_v = A[*v_ptr].get_position();
if (pos_u == pos_v)
{//if2 (Exception handling if two nodes have the same position)
pos_u = N.choose_distinct_random_point_in_radius_epsilon(pos_u);
}//if2
vector_v_minus_u = pos_v - pos_u;
norm_v_minus_u = vector_v_minus_u.norm();
if(!N.f_rep_near_machine_precision(norm_v_minus_u,f_rep_u_on_v))
{
scalar = f_rep_scalar(norm_v_minus_u)/norm_v_minus_u ;
f_rep_u_on_v.m_x = scalar * vector_v_minus_u.m_x;
f_rep_u_on_v.m_y = scalar * vector_v_minus_u.m_y;
}
F_direct[*v_ptr] = F_direct[*v_ptr] + f_rep_u_on_v;
F_direct[*u_ptr] = F_direct[*u_ptr] - f_rep_u_on_v;
}//for
}//if
}//forall2
//Step 3: calculated forces to nodes in act_contained_nodes() of
//leaf_ptr->get_D2()
(*act_leaf_ptr)->get_D2(non_neighboured_leaves);
forall_listiterators(QuadTreeNodeNM*, non_neighbour_leaf_ptr,
non_neighboured_leaves)
{//forall3
(*non_neighbour_leaf_ptr)->get_contained_nodes(
non_neighbour_contained_nodes);
forall_listiterators(node,v_ptr,act_contained_nodes)
forall_listiterators(node, u_ptr,non_neighbour_contained_nodes)
{//for
pos_u = A[*u_ptr].get_position();
pos_v = A[*v_ptr].get_position();
if (pos_u == pos_v)
{//if2 (Exception handling if two nodes have the same position)
pos_u = N.choose_distinct_random_point_in_radius_epsilon(pos_u);
}//if2
vector_v_minus_u = pos_v - pos_u;
norm_v_minus_u = vector_v_minus_u.norm();
if(!N.f_rep_near_machine_precision(norm_v_minus_u,f_rep_u_on_v))
{
scalar = f_rep_scalar(norm_v_minus_u)/norm_v_minus_u ;
f_rep_u_on_v.m_x = scalar * vector_v_minus_u.m_x;
f_rep_u_on_v.m_y = scalar * vector_v_minus_u.m_y;
}
F_direct[*v_ptr] = F_direct[*v_ptr] + f_rep_u_on_v;
}//for
}//forall3
}//if(usual case)
else //special case (more then particles_in_leaves() particles in this leaf)
{//else
forall_listiterators(node, v_ptr, act_contained_nodes)
{
pos_v = A[*v_ptr].get_position();
pos_u = N.choose_distinct_random_point_in_radius_epsilon(pos_v);
vector_v_minus_u = pos_v - pos_u;
norm_v_minus_u = vector_v_minus_u.norm();
if(!N.f_rep_near_machine_precision(norm_v_minus_u,f_rep_u_on_v))
{
scalar = f_rep_scalar(norm_v_minus_u)/norm_v_minus_u ;
f_rep_u_on_v.m_x = scalar * vector_v_minus_u.m_x;
f_rep_u_on_v.m_y = scalar * vector_v_minus_u.m_y;
}
F_direct[*v_ptr] = F_direct[*v_ptr] + f_rep_u_on_v;
}
}//else
}//forall
}
inline void NMM::add_rep_forces(
const Graph& G,
NodeArray& F_direct,
NodeArray& F_multipole_exp,
NodeArray& F_local_exp,
NodeArray& F_rep)
{
node v;
forall_nodes(v,G)
{
F_rep[v] = F_direct[v]+F_local_exp[v]+F_multipole_exp[v];
}
}
inline double NMM::f_rep_scalar(double d)
{
if (d > 0)
{
return 1/d;
}
else
{
cout<<"Error NMM:: f_rep_scalar nodes at same position"<