http://e6.ijs.si/medusa/wiki/api.php?action=feedcontributions&user=JLapajne&feedformat=atomMedusa: Coordinate Free Mehless Method implementation - User contributions [en]2026-05-08T16:32:25ZUser contributionsMediaWiki 1.27.1http://e6.ijs.si/medusa/wiki/index.php?title=Analysis_of_MLSM_performance&diff=1139Analysis of MLSM performance2017-05-04T18:56:58Z<p>JLapajne: /* Convergence with respect to number of nodes */</p>
<hr />
<div>=Solving Diffusion equation=<br />
For starters, we can solve simple Diffusion equation<br />
$ \nabla^2 u = \frac{\partial u}{\partial t} $.<br />
<br />
We solved the equation on a square $\Omega = [0, a] \times [0, a]$ with<br />
Dirichlet boundary conditions $ \left. u\right|_{\partial \Omega} = 0 $ and<br />
initial state $ u(t = 0) = 1$.<br />
<br />
An analytical solution for this domain is known, and we use it to evaluate or<br />
own solution. <br />
\begin{equation} <br />
u(\vec{p}, t) = \sum_{\substack{n=1 \\ n \text{<br />
odd}}}^\infty\sum_{\substack{m=1 \\ m \text{ odd}}}^\infty \frac{1}{\pi^2}<br />
\frac{16 a^2}{nm} \sin\left(\frac{\pi n}{a}p_x\right) \sin\left(\frac{\pi<br />
m}{a}p_y\right) e^{-\frac{\pi^2 (n^2+m^2)}{a^2}t} <br />
\end{equation}<br />
Because the solution is<br />
given in the series form, we only compare to the finite approximation, summing<br />
to $N = 100$ instead of infinity. Solution is on <xr id="fig:square_heat"/>.<br />
See the code for solving diffusion [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp here].<br />
<br />
<figure id="fig:square_heat"><br />
[[File:square_heat.png|thumb|upright=2|center|alt=A square of nodes coloured according to the solution(with smaller and larger node density)|<caption>A picture of our solution (with smaller and larger node density)</caption>]]<br />
</figure><br />
<br />
==Convergence with respect to number of nodes==<br />
<br />
<figure id="fig:node_convergence"><br />
[[File:node_convergence.png|thumb|upright=2|center|alt=Graph of errors with respect to number of nodes|<caption>Convergence with respect to number of nodes</caption>]]<br />
</figure><br />
<br />
We tested the method with a fixed time step of $ \Delta t = 1\cdot 10^{-5}$<br />
on a unit square ($a = 1$). Results are on <xr id="fig:node_convergence"/>. Monomial basis of $6$ monomials was used and $12$<br />
closest nodes counted as support for each node. After more than $250$ nodes of<br />
discretization in each dimension the method diverges, which is expected. The<br />
stability criterion for diffusion equation in two dimensions is $\Delta t \leq<br />
\frac{1}{4} \Delta x^2$, where $\Delta x$ is the spatial discretization<br />
step in one dimension. In our case, at $250$ nodes per side, the right hand side<br />
yields $\frac{1}{4}\cdot\frac{1}{250}\cdot\frac{1}{250} = 4\times 10^{-6}$,<br />
so our method is stable within the expected region.<br />
<br />
<figure id="fig:node_convergence_5"><br />
[[File:node_convergence_5.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for Gaussian and Monomial basis</caption>]]<br />
</figure><br />
<br />
On <xr id="fig:node_convergence_5"/> is another image of convergence, this time using monomial basis $\{1, x,<br />
y, x^2, y^2\}$ and Gaussian basis with discretization step $ \Delta t = 5<br />
\cdot 10^{-6}$ and 5 support nodes. Total node count was $N = 2500$.<br />
Error was calculated after $0.01$ time units have elapsed.<br />
<br />
<figure id="fig:node_convergence_comparison"><br />
[[File:explicit_implicit_diffusion_comparison_ver2.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for explicit and implicit mixed boundary condition (Neumann is enforced in two different ways).</caption>]]<br />
</figure><br />
<br />
Figure 4 shows convergence of 1/4 of the whole solution with Dirichlet boundary conditions (Fig. 1). The nodes were spanning the area $[0, a/2]^2$. As expected, higher number of nodes yields lower max error up to some maximal number of nodes. The latter is lower for explicit than for implicit method.<br />
<br />
==Convergence with respect to number of time steps==<br />
<br />
<figure id="fig:timestep_convergence"><br />
[[File:timestep_convergence.png|thumb|upright=2|center|<caption>Convergence with respect to different timesteps</caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count with different time steps on the same<br />
domain as above. Results are on <xr id="fig:timestep_convergence"/>. For large time steps the method diverges, but once it starts<br />
converging the precision increases steadily as the time step decreases, until it<br />
hits its lower limit. This behaviour is expected. The error was calculated<br />
against the analytical solution above after $0.005$ units of time have passed. A<br />
monomial basis up to order $2$ inclusive ($m = 6$) was used and the support<br />
size was $n = 12$.<br />
<br />
==Using Gaussian basis==<br />
<br />
<figure id="fig:gauss_sigma_dependence"><br />
[[File:gauss_sigma_dependence.png|thumb|upright=2|center|<caption> Graph of error with respect to Gaussian basis $\sigma$ </caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count of $N = 2500$ with spatial step<br />
discretization of $\Delta x = \frac{1}{50}$. We used the Gaussian basis of<br />
$m = 5$ functions with support size $n = 13$. Error was calculated<br />
against analytical solution above after $0.01$ time units. A time step of<br />
$\Delta t = 10^{-5}$ was used. Results are on <xr id="fig:gauss_sigma_dependence"/><br />
<br />
As we can see from the graph, there exists an interval where the choice of<br />
$\sigma$ does not matter much, but outside of this interval, the method diverges<br />
rapidly. Care from the user side must be taken, to choose $\sigma$<br />
appropriately, with respect to plot above.<br />
<br />
== Convergence with respect to shape parameters ==<br />
<br />
Choice of shape of the basis function can have significant effect on solution of the system.<br />
We tested the space of feasible shape parameters on the diffusion equation with $\Delta x = \Delta y = 0.01$<br />
and $m = n = 9$. The weight function does not have much influence here, except on the condition <br />
number and any singular values that get cut off because of that.<br />
<br />
Results for Gaussian, MQ and IMQ basis functions are presented in <xr id="fig:shape_space_gau"/>, <xr id="fig:shape_space_mq"/> and <xr id="fig:shape_space_imq"/>.<br />
<br />
<figure id="fig:shape_space_gau"><br />
[[File:sigma_depedance_error_gau.png|thumb|left|800px|<caption>Error as a function of Gaussian basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_mq"><br />
[[File:sigma_depedance_error_imq.png|thumb|800px|<caption>Error as a function of MQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mq.png|thumb|upright=1|left|800px|<caption>Error as a function of IMQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mon.png|thumb|upright=1|800px|<caption>Error as a function of Gaussian weight's $\sigma_W$ for Monomials of order 2.</caption>]]<br />
</figure><br />
<br />
<div style="clear:both"></div><br />
<br />
==Solving in 3D==<br />
<br />
A 3-dimensional case on domain $[0, 1]^3$ was tested on $N = 20^3$<br />
nodes, making the discretization step $\Delta x = 0.05$. Support size of<br />
$n=10$ with $m=10$ Gaussian basis functions was used. Their<br />
normalization parameter was $\sigma = 60\Delta x$. A time step of $\Delta<br />
t = 10^{-5}$ and an explicit Euler method was used to calculate the solution<br />
of to $0.01$ time units. Resulting function is on <xr id="fig:diffusion3d"/>.<br />
<br />
<br />
<figure id="fig:diffusion3d"><br />
[[File:diffusion3d.png|thumb|upright =2|center|<caption>A $ 3 $-dimensional soulution with an explicit Euler method</caption>]]<br />
</figure><br />
<br />
==Solving Dirichlet with [http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classEngineMLS.html EngineMLS]==<br />
Example code using explicit stepping to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp example code]<br />
<br />
==Solving Dirichlet with MLSM operators==<br />
Example code using explicit stepping and MLSM operators to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion_mlsm_operators.cpp example code]<br />
<br />
==Solving mixed boundary conditions with MLSM operators==<br />
By using MLSM operators and utilizing its `[http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classmm_1_1MLSM.html MLSM::neumann]` method we can also solve<br />
partial diffusion equations. Solution is on <xr id="fig:quarter_diffusion"/>.<br />
<br />
<figure id="fig:quarter_diffusion"><br />
[[File:quarter_diffusion.png|400px|thumb|upright=2|center|<caption>Example of solving using Neumann boundary conditions </caption>]]<br />
</figure><br />
<br />
Example code showing the use of Neumann boundary conditions: [https://gitlab.com/e62Lab/e62numcodes/blob/master/test/mlsm_operators_test.cpp example code]<br />
<br />
We also support more -- interesting -- domains :) On <xr id="fig:mikimouse_heat"/> we see a domain in a shape of Miki mouse.<br />
<br />
<figure id="fig:mikimouse_heat"><br />
[[File:mikimouse_heat.png|400px|thumb|upright=2|center|<caption>Miki mouse domain </caption>]]<br />
</figure><br />
<br />
<!--<br />
==Convergence analysis for one dimensional diffusion equation==<br />
We now tried to solve the diffusion equation in one dimension <br />
$u_t = u_{xx}$ on the interval $[0,1]$<br />
with Dirichlet boundary conditions<br />
$u(0, t) = u(1,t) = 0 $<br />
and an initial state<br />
$u(x,0) = \sin(\pi x).$<br />
<br />
An analytical solution for this domain is given in closed form <br />
\begin{equation}<br />
u(x,t) = \sin(\pi x) \exp(-\pi^2 t),<br />
\end{equation}<br />
so we can use it to evaluate our approximation.<br />
<br />
We tested the method with a fixed time step of $\Delta t = 1\cdot 10^{-8}$ on a unit interval, that was discretized into $ N $ nodes. Monomial basis of $m$ monomials was used and $n$ closest nodes counted as support for each node. <br />
<br />
<strong> Note: </strong> the results are <strong>not confirmed </strong> and serve merely as an orientation. Errors may have occured. <br />
<br />
==Convergence with respect to number of nodes and support size==<br />
<br />
Convergence of the method was tested on different number of nodes $N = 15, 20, 25, \ldots, 200$. First we tested the method, when support size is equal to number of monomials in the basis ($n =m$), which we call local interpolation, for different values of $n = 3, 5, 7$. We then compared that to local approximation, where we take $n = 12$ supporting nodes and we only change number of monomials. We used gaussian weight function with $\sigma = 1\cdot \Delta x$. Error was calculated after $0.001$ time units have elapsed. We can see the results on <xr id="fig:diffusionConvergence1d"/>.<br />
<br />
<figure id="fig:diffusionConvergence1d"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional ddiffusion equation using monomial basis.</caption>]]<br />
</figure><br />
<br />
Larger is the support size, faster is convergence. This rule applies for local interpolation as well as approximation. But when number of nodes is too large, the method diverges faster for larger support sizes. Error for approximation is bigger than for interpolation (if we compare it for same basis size), but the method diverges with interpolation for smaller number of nodes.<br />
<br />
==Convergence with respect to weight function==<br />
<br />
We also checked the convergence for different weight functions. We tested it on number of nodes $N = 15, 20, 25, \ldots, 200$, with monomial basis with $m = 5$, support size $n = 12$ and gaussian weight function with $\sigma = w\cdot \Delta x$ for $w = 0.1$. Error was calculated after $0.001$ time units have elapsed with $\Delta t = 10^{-8}$. We can see the results on <xr id="fig:diffusionConvergence1dWeights"/>.<br />
<br />
<figure id="fig:diffusionConvergence1dWeights"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12_weights.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional diffusion equation using monomial basis with respect to weight function. </caption>]]<br />
</figure><br />
<br />
When the $\sigma$ is too small, other nodes on the support domain are not taken into consideration, so taking $n = 12$ or $n = 3$ would give the same results. However when $\sigma$ is too large, all the nodes on the support domain are contributing and there is no significant difference between different weight functions.<br />
--><br />
<br />
= Solving Electrostatics =<br />
<br />
This example is taken from the [http://www.freefem.org/ff++/ftp/freefem++doc.pdf FreeFem++ manual (page 235)].<br />
<br />
Assuming there is no current and the charge distribution is time independent, the electric field $\b{E}$ satisfies<br />
\begin{equation}\label{eq:electrostatics}<br />
\b{\nabla}\cdot\b{E} = \frac{\rho}{\epsilon}, \quad \b{\nabla} \times \b{E} = 0<br />
\end{equation}<br />
where $\rho$ is the charge density and $\epsilon$ is the permittivity. If we introduce an electrostatics potential $\phi$ such that<br />
\begin{equation}<br />
\b{E} = -\b{\nabla}\phi,<br />
\end{equation}<br />
we can insert it into the first equation in (\ref{eq:electrostatics}) resulting in Poisson's equation<br />
\begin{equation}\label{eq:poisson}<br />
\b{\nabla}^2 \phi = -\frac{\rho}{\epsilon}.<br />
\end{equation}<br />
In the absence of unpaired electric charge equation (\ref{eq:poisson}) becomes Laplace's equation<br />
\begin{equation}<br />
\b{\nabla}^2 \phi = 0<br />
\end{equation}<br />
<br />
We now solve this equation for a circular enclosure with two rectangular holes. The boundary of the circular enclosure is held at constant potential $0$ V. The two rectangular holes are held at constant potentials $+1$ V and $-1$ V, respectively. A coloured scatter plot is available in figure <xr id="fig:electro_statics"/>.<br />
<br />
<br />
<div class="toccolours mw-collapsible mw-collapsed"> '''Code for electrostatics problem'''<br />
<div class="mw-collapsible-content"><br />
<br />
<syntaxhighlight lang="c++" line><br />
// domain parameters, size, support and monomial order<br />
int domain_size = 3000;<br />
size_t n = 15; <br />
size_t m = 3; <br />
<br />
double radius = 5.0;<br />
double width = 0.3;<br />
double height = 3.0;<br />
<br />
// extra boundary labels<br />
int LEFT = -10;<br />
int RIGHT = -20;<br />
<br />
// build circular domain<br />
CircleDomain<Vec2d> domain({0,0},radius);<br />
domain.fillUniformInterior(domain_size);<br />
domain.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
<br />
// build left rectangle<br />
RectangleDomain<Vec2d> left({-2-width,-height},{-2+width,height});<br />
left.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
left.types[left.types == BOUNDARY] = LEFT;<br />
domain.subtract(left);<br />
<br />
// build right rectangle<br />
RectangleDomain<Vec2d> right({2-width,-height},{2+width,height});<br />
right.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
right.types[right.types == BOUNDARY] = RIGHT;<br />
domain.subtract(right);<br />
<br />
// relax domain and find supports<br />
domain.relax(50, 1e-2, 1.0, 3, 1000);<br />
domain.findSupport(n);<br />
<br />
// get interior and boundary ranges<br />
Range<int> interior = domain.types == INTERNAL;<br />
Range<int> boundary = domain.types == BOUNDARY;<br />
Range<int> leftrect = domain.types == LEFT;<br />
Range<int> rightrect = domain.types == RIGHT;<br />
<br />
// initialize unknown concentration field and right-hand side<br />
VecXd phi(domain.size(),1);<br />
VecXd RHS(domain.size(),1);<br />
<br />
// initialize interior values (this is an important step)<br />
RHS[interior] = 0.0;<br />
<br />
// set dirichlet boundary conditions<br />
RHS[boundary] = 0.0;<br />
RHS[leftrect] = -1.0;<br />
RHS[rightrect] = 1.0;<br />
<br />
// prepare shape functions and laplacian<br />
std::vector<Triplet<double>> l_coeff;<br />
for (auto& c : interior) {<br />
Range<Vec2d> supp_domain = domain.positions[domain.support[c]];<br />
EngineMLS<Vec2d, Monomials, Gaussians> MLS(m,supp_domain,pow(domain.characteristicDistance(),2));<br />
VecXd shape = MLS.getShapeAt(supp_domain[0], {{2, 0}}) +<br />
MLS.getShapeAt(supp_domain[0], {{0, 2}});<br />
for (size_t i = 0; i < supp_domain.size(); ++i) {<br />
l_coeff.emplace_back(c, domain.support[c][i], shape[i]);<br />
}<br />
}<br />
<br />
// prepare dirichlet boundaries<br />
std::vector<Triplet<double>> b_coeff;<br />
for (auto& c : domain.types < 0) {<br />
b_coeff.emplace_back(c, c, 1.0);<br />
}<br />
<br />
// prepare matrices<br />
SparseMatrix<double> laplace_m(domain.size(), domain.size());<br />
laplace_m.setFromTriplets(l_coeff.begin(), l_coeff.end());<br />
SparseMatrix<double> boundary_m(domain.size(), domain.size());<br />
boundary_m.setFromTriplets(b_coeff.begin(), b_coeff.end());<br />
<br />
// draw<br />
std::thread th([&] { draw2D(domain.positions, phi); });<br />
<br />
// initialize solver<br />
Eigen::BiCGSTAB<SparseMatrix<double>> solver;<br />
<br />
// join matrices together<br />
SparseMatrix<double> tmp = boundary_m + laplace_m;<br />
solver.compute(tmp);<br />
<br />
// solve system of equations<br />
phi = solver.solve(RHS);<br />
std::cout << "#iterations: " << solver.iterations() << std::endl;<br />
std::cout << "estimated error: " << solver.error() << std::endl;<br />
<br />
// end drawing<br />
th.join(); <br />
</syntaxhighlight><br />
<br />
</div><br />
</div><br />
<br />
<br />
<figure id="fig:electro_statics"><br />
[[File:electrostatics.png|thumb|upright=2|center|<caption> Simple electrostatics problem solved by implicit method with 15 node monomial supports and gaussian weight functions. </caption>]]<br />
</figure><br />
<br />
= Convection-diffusion =<br />
<br />
The following example problem is adapted from [http://servidor.demec.ufpr.br/CFD/bibliografia/Ferziger%20Peric%20-%20Computational%20Methods%20for%20Fluid%20Dynamics,%203rd%20Ed%20-%202002.pdf Ferziger & Perič (2002) (pages 82-86)].<br />
<br />
<figure id="fig:fvm_problem"><br />
[[File:Screenshot_2016-11-28_11-05-03.png|thumb|400px|upright=2|centre|<caption> Geometry and boundary conditions for the scalar transport in a stagnation point flow. </caption>]]<br />
</figure><br />
<br />
We now consider the problem of a scalar quantity transported in a known velocity field $\b{u}$. The velocity field is given by<br />
\[\b{u} = (u_x,u_y) = (x, -y),\]<br />
which represents the flow near a stagnation point $(x,y) = (0,0)$. The streamlines $xy=\mathrm{const.}$ and change direction with respect to the Cartesian grid. The convection-diffusion equation (also known as the scalar transport equation) to be solved is<br />
\[\alpha \b{\nabla}^2 \phi - \b{u}\cdot \b{\nabla}\phi = 0\]<br />
with the following boundary conditions:<br />
* $\phi = 0$ along the north (inlet) boundary;<br />
* linear variation $\phi = (1 - y)$ along the west boundary ($x = 0$);<br />
* symmetry condition on the south boundary;<br />
* zero gradient at the east (outlet) boundary.<br />
The parameter $\alpha$ is the diffusivity for the scalar quantity $\phi$. The problem and boundary conditions are schematically represented in the <xr id="fig:fvm_problem"/>. The solutions obtained by Ferziger and Perić are shown in <xr id="fig:fvm_problem_solution"/>.<br />
<br />
<figure id="fig:fvm_problem_solution"><br />
[[File:Screenshot_2016-11-28_11-07-44.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right). </caption>]]<br />
</figure><br />
<br />
The meshless solution is obtained on regular point arrangements of sizes 100 × 100 and 300 × 300, for the 2 cases of diffusivity $\alpha = 0.01$ and $\alpha = 0.001$, respectively. For comparison we have plotted the contour lines at the same intervals as the original authors of this example (see <xr id="fig:fvm_problem_meshless"/>). The basis functions are the 5 monomials $1$, $x$, $y$, $x^2$, $y^2$.<br />
<br />
<figure id="fig:fvm_problem_meshless"><br />
[[File:fvm_convection.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right), obtaines with meshless solver. </caption>]]<br />
</figure></div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=File:Explicit_implicit_diffusion_comparison_ver2.png&diff=1138File:Explicit implicit diffusion comparison ver2.png2017-05-04T18:56:07Z<p>JLapajne: File uploaded with MsUpload</p>
<hr />
<div>File uploaded with MsUpload</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=Analysis_of_MLSM_performance&diff=1127Analysis of MLSM performance2017-04-24T09:02:18Z<p>JLapajne: /* Convergence with respect to number of nodes */</p>
<hr />
<div>=Solving Diffusion equation=<br />
For starters, we can solve simple Diffusion equation<br />
$ \nabla^2 u = \frac{\partial u}{\partial t} $.<br />
<br />
We solved the equation on a square $\Omega = [0, a] \times [0, a]$ with<br />
Dirichlet boundary conditions $ \left. u\right|_{\partial \Omega} = 0 $ and<br />
initial state $ u(t = 0) = 1$.<br />
<br />
An analytical solution for this domain is known, and we use it to evaluate or<br />
own solution. <br />
\begin{equation} <br />
u(\vec{p}, t) = \sum_{\substack{n=1 \\ n \text{<br />
odd}}}^\infty\sum_{\substack{m=1 \\ m \text{ odd}}}^\infty \frac{1}{\pi^2}<br />
\frac{16 a^2}{nm} \sin\left(\frac{\pi n}{a}p_x\right) \sin\left(\frac{\pi<br />
m}{a}p_y\right) e^{-\frac{\pi^2 (n^2+m^2)}{a^2}t} <br />
\end{equation}<br />
Because the solution is<br />
given in the series form, we only compare to the finite approximation, summing<br />
to $N = 100$ instead of infinity. Solution is on <xr id="fig:square_heat"/>.<br />
See the code for solving diffusion [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp here].<br />
<br />
<figure id="fig:square_heat"><br />
[[File:square_heat.png|thumb|upright=2|center|alt=A square of nodes coloured according to the solution(with smaller and larger node density)|<caption>A picture of our solution (with smaller and larger node density)</caption>]]<br />
</figure><br />
<br />
==Convergence with respect to number of nodes==<br />
<br />
<figure id="fig:node_convergence"><br />
[[File:node_convergence.png|thumb|upright=2|center|alt=Graph of errors with respect to number of nodes|<caption>Convergence with respect to number of nodes</caption>]]<br />
</figure><br />
<br />
We tested the method with a fixed time step of $ \Delta t = 1\cdot 10^{-5}$<br />
on a unit square ($a = 1$). Results are on <xr id="fig:node_convergence"/>. Monomial basis of $6$ monomials was used and $12$<br />
closest nodes counted as support for each node. After more than $250$ nodes of<br />
discretization in each dimension the method diverges, which is expected. The<br />
stability criterion for diffusion equation in two dimensions is $\Delta t \leq<br />
\frac{1}{4} \Delta x^2$, where $\Delta x$ is the spatial discretization<br />
step in one dimension. In our case, at $250$ nodes per side, the right hand side<br />
yields $\frac{1}{4}\cdot\frac{1}{250}\cdot\frac{1}{250} = 4\times 10^{-6}$,<br />
so our method is stable within the expected region.<br />
<br />
<figure id="fig:node_convergence_5"><br />
[[File:node_convergence_5.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for Gaussian and Monomial basis</caption>]]<br />
</figure><br />
<br />
On <xr id="fig:node_convergence_5"/> is another image of convergence, this time using monomial basis $\{1, x,<br />
y, x^2, y^2\}$ and Gaussian basis with discretization step $ \Delta t = 5<br />
\cdot 10^{-6}$ and 5 support nodes. Total node count was $N = 2500$.<br />
Error was calculated after $0.01$ time units have elapsed.<br />
<br />
<figure id="fig:node_convergence_comparison"><br />
[[File:explicit_implicit_diffusion_comparison.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for explicit and implicit mixed boundary condition (Neumann is enforced in two different ways).</caption>]]<br />
</figure><br />
<br />
Figure 4 shows convergence of 1/4 of the whole solution with Dirichlet boundary conditions (Fig. 1). The nodes were spanning the area $[0, 0.5]^2$. As expected, higher number of nodes yields lower max error up to some maximal number of nodes. The latter is lower for explicit than for implicit method.<br />
<br />
==Convergence with respect to number of time steps==<br />
<br />
<figure id="fig:timestep_convergence"><br />
[[File:timestep_convergence.png|thumb|upright=2|center|<caption>Convergence with respect to different timesteps</caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count with different time steps on the same<br />
domain as above. Results are on <xr id="fig:timestep_convergence"/>. For large time steps the method diverges, but once it starts<br />
converging the precision increases steadily as the time step decreases, until it<br />
hits its lower limit. This behaviour is expected. The error was calculated<br />
against the analytical solution above after $0.005$ units of time have passed. A<br />
monomial basis up to order $2$ inclusive ($m = 6$) was used and the support<br />
size was $n = 12$.<br />
<br />
==Using Gaussian basis==<br />
<br />
<figure id="fig:gauss_sigma_dependence"><br />
[[File:gauss_sigma_dependence.png|thumb|upright=2|center|<caption> Graph of error with respect to Gaussian basis $\sigma$ </caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count of $N = 2500$ with spatial step<br />
discretization of $\Delta x = \frac{1}{50}$. We used the Gaussian basis of<br />
$m = 5$ functions with support size $n = 13$. Error was calculated<br />
against analytical solution above after $0.01$ time units. A time step of<br />
$\Delta t = 10^{-5}$ was used. Results are on <xr id="fig:gauss_sigma_dependence"/><br />
<br />
As we can see from the graph, there exists an interval where the choice of<br />
$\sigma$ does not matter much, but outside of this interval, the method diverges<br />
rapidly. Care from the user side must be taken, to choose $\sigma$<br />
appropriately, with respect to plot above.<br />
<br />
== Convergence with respect to shape parameters ==<br />
<br />
Choice of shape of the basis function can have significant effect on solution of the system.<br />
We tested the space of feasible shape parameters on the diffusion equation with $\Delta x = \Delta y = 0.01$<br />
and $m = n = 9$. The weight function does not have much influence here, except on the condition <br />
number and any singular values that get cut off because of that.<br />
<br />
Results for Gaussian, MQ and IMQ basis functions are presented in <xr id="fig:shape_space_gau"/>, <xr id="fig:shape_space_mq"/> and <xr id="fig:shape_space_imq"/>.<br />
<br />
<figure id="fig:shape_space_gau"><br />
[[File:sigma_depedance_error_gau.png|thumb|left|800px|<caption>Error as a function of Gaussian basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_mq"><br />
[[File:sigma_depedance_error_imq.png|thumb|800px|<caption>Error as a function of MQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mq.png|thumb|upright=1|left|800px|<caption>Error as a function of IMQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mon.png|thumb|upright=1|800px|<caption>Error as a function of Gaussian weight's $\sigma_W$ for Monomials of order 2.</caption>]]<br />
</figure><br />
<br />
<div style="clear:both"></div><br />
<br />
==Solving in 3D==<br />
<br />
A 3-dimensional case on domain $[0, 1]^3$ was tested on $N = 20^3$<br />
nodes, making the discretization step $\Delta x = 0.05$. Support size of<br />
$n=10$ with $m=10$ Gaussian basis functions was used. Their<br />
normalization parameter was $\sigma = 60\Delta x$. A time step of $\Delta<br />
t = 10^{-5}$ and an explicit Euler method was used to calculate the solution<br />
of to $0.01$ time units. Resulting function is on <xr id="fig:diffusion3d"/>.<br />
<br />
<br />
<figure id="fig:diffusion3d"><br />
[[File:diffusion3d.png|thumb|upright =2|center|<caption>A $ 3 $-dimensional soulution with an explicit Euler method</caption>]]<br />
</figure><br />
<br />
==Solving Dirichlet with [http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classEngineMLS.html EngineMLS]==<br />
Example code using explicit stepping to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp example code]<br />
<br />
==Solving Dirichlet with MLSM operators==<br />
Example code using explicit stepping and MLSM operators to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion_mlsm_operators.cpp example code]<br />
<br />
==Solving mixed boundary conditions with MLSM operators==<br />
By using MLSM operators and utilizing its `[http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classmm_1_1MLSM.html MLSM::neumann]` method we can also solve<br />
partial diffusion equations. Solution is on <xr id="fig:quarter_diffusion"/>.<br />
<br />
<figure id="fig:quarter_diffusion"><br />
[[File:quarter_diffusion.png|400px|thumb|upright=2|center|<caption>Example of solving using Neumann boundary conditions </caption>]]<br />
</figure><br />
<br />
Example code showing the use of Neumann boundary conditions: [https://gitlab.com/e62Lab/e62numcodes/blob/master/test/mlsm_operators_test.cpp example code]<br />
<br />
We also support more -- interesting -- domains :) On <xr id="fig:mikimouse_heat"/> we see a domain in a shape of Miki mouse.<br />
<br />
<figure id="fig:mikimouse_heat"><br />
[[File:mikimouse_heat.png|400px|thumb|upright=2|center|<caption>Miki mouse domain </caption>]]<br />
</figure><br />
<br />
<!--<br />
==Convergence analysis for one dimensional diffusion equation==<br />
We now tried to solve the diffusion equation in one dimension <br />
$u_t = u_{xx}$ on the interval $[0,1]$<br />
with Dirichlet boundary conditions<br />
$u(0, t) = u(1,t) = 0 $<br />
and an initial state<br />
$u(x,0) = \sin(\pi x).$<br />
<br />
An analytical solution for this domain is given in closed form <br />
\begin{equation}<br />
u(x,t) = \sin(\pi x) \exp(-\pi^2 t),<br />
\end{equation}<br />
so we can use it to evaluate our approximation.<br />
<br />
We tested the method with a fixed time step of $\Delta t = 1\cdot 10^{-8}$ on a unit interval, that was discretized into $ N $ nodes. Monomial basis of $m$ monomials was used and $n$ closest nodes counted as support for each node. <br />
<br />
<strong> Note: </strong> the results are <strong>not confirmed </strong> and serve merely as an orientation. Errors may have occured. <br />
<br />
==Convergence with respect to number of nodes and support size==<br />
<br />
Convergence of the method was tested on different number of nodes $N = 15, 20, 25, \ldots, 200$. First we tested the method, when support size is equal to number of monomials in the basis ($n =m$), which we call local interpolation, for different values of $n = 3, 5, 7$. We then compared that to local approximation, where we take $n = 12$ supporting nodes and we only change number of monomials. We used gaussian weight function with $\sigma = 1\cdot \Delta x$. Error was calculated after $0.001$ time units have elapsed. We can see the results on <xr id="fig:diffusionConvergence1d"/>.<br />
<br />
<figure id="fig:diffusionConvergence1d"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional ddiffusion equation using monomial basis.</caption>]]<br />
</figure><br />
<br />
Larger is the support size, faster is convergence. This rule applies for local interpolation as well as approximation. But when number of nodes is too large, the method diverges faster for larger support sizes. Error for approximation is bigger than for interpolation (if we compare it for same basis size), but the method diverges with interpolation for smaller number of nodes.<br />
<br />
==Convergence with respect to weight function==<br />
<br />
We also checked the convergence for different weight functions. We tested it on number of nodes $N = 15, 20, 25, \ldots, 200$, with monomial basis with $m = 5$, support size $n = 12$ and gaussian weight function with $\sigma = w\cdot \Delta x$ for $w = 0.1$. Error was calculated after $0.001$ time units have elapsed with $\Delta t = 10^{-8}$. We can see the results on <xr id="fig:diffusionConvergence1dWeights"/>.<br />
<br />
<figure id="fig:diffusionConvergence1dWeights"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12_weights.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional diffusion equation using monomial basis with respect to weight function. </caption>]]<br />
</figure><br />
<br />
When the $\sigma$ is too small, other nodes on the support domain are not taken into consideration, so taking $n = 12$ or $n = 3$ would give the same results. However when $\sigma$ is too large, all the nodes on the support domain are contributing and there is no significant difference between different weight functions.<br />
--><br />
<br />
= Solving Electrostatics =<br />
<br />
This example is taken from the [http://www.freefem.org/ff++/ftp/freefem++doc.pdf FreeFem++ manual (page 235)].<br />
<br />
Assuming there is no current and the charge distribution is time independent, the electric field $\b{E}$ satisfies<br />
\begin{equation}\label{eq:electrostatics}<br />
\b{\nabla}\cdot\b{E} = \frac{\rho}{\epsilon}, \quad \b{\nabla} \times \b{E} = 0<br />
\end{equation}<br />
where $\rho$ is the charge density and $\epsilon$ is the permittivity. If we introduce an electrostatics potential $\phi$ such that<br />
\begin{equation}<br />
\b{E} = -\b{\nabla}\phi,<br />
\end{equation}<br />
we can insert it into the first equation in (\ref{eq:electrostatics}) resulting in Poisson's equation<br />
\begin{equation}\label{eq:poisson}<br />
\b{\nabla}^2 \phi = -\frac{\rho}{\epsilon}.<br />
\end{equation}<br />
In the absence of unpaired electric charge equation (\ref{eq:poisson}) becomes Laplace's equation<br />
\begin{equation}<br />
\b{\nabla}^2 \phi = 0<br />
\end{equation}<br />
<br />
We now solve this equation for a circular enclosure with two rectangular holes. The boundary of the circular enclosure is held at constant potential $0$ V. The two rectangular holes are held at constant potentials $+1$ V and $-1$ V, respectively. A coloured scatter plot is available in figure <xr id="fig:electro_statics"/>.<br />
<br />
<br />
<div class="toccolours mw-collapsible mw-collapsed"> '''Code for electrostatics problem'''<br />
<div class="mw-collapsible-content"><br />
<br />
<syntaxhighlight lang="c++" line><br />
// domain parameters, size, support and monomial order<br />
int domain_size = 3000;<br />
size_t n = 15; <br />
size_t m = 3; <br />
<br />
double radius = 5.0;<br />
double width = 0.3;<br />
double height = 3.0;<br />
<br />
// extra boundary labels<br />
int LEFT = -10;<br />
int RIGHT = -20;<br />
<br />
// build circular domain<br />
CircleDomain<Vec2d> domain({0,0},radius);<br />
domain.fillUniformInterior(domain_size);<br />
domain.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
<br />
// build left rectangle<br />
RectangleDomain<Vec2d> left({-2-width,-height},{-2+width,height});<br />
left.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
left.types[left.types == BOUNDARY] = LEFT;<br />
domain.subtract(left);<br />
<br />
// build right rectangle<br />
RectangleDomain<Vec2d> right({2-width,-height},{2+width,height});<br />
right.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
right.types[right.types == BOUNDARY] = RIGHT;<br />
domain.subtract(right);<br />
<br />
// relax domain and find supports<br />
domain.relax(50, 1e-2, 1.0, 3, 1000);<br />
domain.findSupport(n);<br />
<br />
// get interior and boundary ranges<br />
Range<int> interior = domain.types == INTERNAL;<br />
Range<int> boundary = domain.types == BOUNDARY;<br />
Range<int> leftrect = domain.types == LEFT;<br />
Range<int> rightrect = domain.types == RIGHT;<br />
<br />
// initialize unknown concentration field and right-hand side<br />
VecXd phi(domain.size(),1);<br />
VecXd RHS(domain.size(),1);<br />
<br />
// initialize interior values (this is an important step)<br />
RHS[interior] = 0.0;<br />
<br />
// set dirichlet boundary conditions<br />
RHS[boundary] = 0.0;<br />
RHS[leftrect] = -1.0;<br />
RHS[rightrect] = 1.0;<br />
<br />
// prepare shape functions and laplacian<br />
std::vector<Triplet<double>> l_coeff;<br />
for (auto& c : interior) {<br />
Range<Vec2d> supp_domain = domain.positions[domain.support[c]];<br />
EngineMLS<Vec2d, Monomials, Gaussians> MLS(m,supp_domain,pow(domain.characteristicDistance(),2));<br />
VecXd shape = MLS.getShapeAt(supp_domain[0], {{2, 0}}) +<br />
MLS.getShapeAt(supp_domain[0], {{0, 2}});<br />
for (size_t i = 0; i < supp_domain.size(); ++i) {<br />
l_coeff.emplace_back(c, domain.support[c][i], shape[i]);<br />
}<br />
}<br />
<br />
// prepare dirichlet boundaries<br />
std::vector<Triplet<double>> b_coeff;<br />
for (auto& c : domain.types < 0) {<br />
b_coeff.emplace_back(c, c, 1.0);<br />
}<br />
<br />
// prepare matrices<br />
SparseMatrix<double> laplace_m(domain.size(), domain.size());<br />
laplace_m.setFromTriplets(l_coeff.begin(), l_coeff.end());<br />
SparseMatrix<double> boundary_m(domain.size(), domain.size());<br />
boundary_m.setFromTriplets(b_coeff.begin(), b_coeff.end());<br />
<br />
// draw<br />
std::thread th([&] { draw2D(domain.positions, phi); });<br />
<br />
// initialize solver<br />
Eigen::BiCGSTAB<SparseMatrix<double>> solver;<br />
<br />
// join matrices together<br />
SparseMatrix<double> tmp = boundary_m + laplace_m;<br />
solver.compute(tmp);<br />
<br />
// solve system of equations<br />
phi = solver.solve(RHS);<br />
std::cout << "#iterations: " << solver.iterations() << std::endl;<br />
std::cout << "estimated error: " << solver.error() << std::endl;<br />
<br />
// end drawing<br />
th.join(); <br />
</syntaxhighlight><br />
<br />
</div><br />
</div><br />
<br />
<br />
<figure id="fig:electro_statics"><br />
[[File:electrostatics.png|thumb|upright=2|center|<caption> Simple electrostatics problem solved by implicit method with 15 node monomial supports and gaussian weight functions. </caption>]]<br />
</figure><br />
<br />
= Convection-diffusion =<br />
<br />
The following example problem is adapted from [http://servidor.demec.ufpr.br/CFD/bibliografia/Ferziger%20Peric%20-%20Computational%20Methods%20for%20Fluid%20Dynamics,%203rd%20Ed%20-%202002.pdf Ferziger & Perič (2002) (pages 82-86)].<br />
<br />
<figure id="fig:fvm_problem"><br />
[[File:Screenshot_2016-11-28_11-05-03.png|thumb|400px|upright=2|centre|<caption> Geometry and boundary conditions for the scalar transport in a stagnation point flow. </caption>]]<br />
</figure><br />
<br />
We now consider the problem of a scalar quantity transported in a known velocity field $\b{u}$. The velocity field is given by<br />
\[\b{u} = (u_x,u_y) = (x, -y),\]<br />
which represents the flow near a stagnation point $(x,y) = (0,0)$. The streamlines $xy=\mathrm{const.}$ and change direction with respect to the Cartesian grid. The convection-diffusion equation (also known as the scalar transport equation) to be solved is<br />
\[\alpha \b{\nabla}^2 \phi - \b{u}\cdot \b{\nabla}\phi = 0\]<br />
with the following boundary conditions:<br />
* $\phi = 0$ along the north (inlet) boundary;<br />
* linear variation $\phi = (1 - y)$ along the west boundary ($x = 0$);<br />
* symmetry condition on the south boundary;<br />
* zero gradient at the east (outlet) boundary.<br />
The parameter $\alpha$ is the diffusivity for the scalar quantity $\phi$. The problem and boundary conditions are schematically represented in the <xr id="fig:fvm_problem"/>. The solutions obtained by Ferziger and Perić are shown in <xr id="fig:fvm_problem_solution"/>.<br />
<br />
<figure id="fig:fvm_problem_solution"><br />
[[File:Screenshot_2016-11-28_11-07-44.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right). </caption>]]<br />
</figure><br />
<br />
The meshless solution is obtained on regular point arrangements of sizes 100 × 100 and 300 × 300, for the 2 cases of diffusivity $\alpha = 0.01$ and $\alpha = 0.001$, respectively. For comparison we have plotted the contour lines at the same intervals as the original authors of this example (see <xr id="fig:fvm_problem_meshless"/>). The basis functions are the 5 monomials $1$, $x$, $y$, $x^2$, $y^2$.<br />
<br />
<figure id="fig:fvm_problem_meshless"><br />
[[File:fvm_convection.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right), obtaines with meshless solver. </caption>]]<br />
</figure></div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=Analysis_of_MLSM_performance&diff=1126Analysis of MLSM performance2017-04-24T09:01:38Z<p>JLapajne: /* Convergence with respect to number of nodes */</p>
<hr />
<div>=Solving Diffusion equation=<br />
For starters, we can solve simple Diffusion equation<br />
$ \nabla^2 u = \frac{\partial u}{\partial t} $.<br />
<br />
We solved the equation on a square $\Omega = [0, a] \times [0, a]$ with<br />
Dirichlet boundary conditions $ \left. u\right|_{\partial \Omega} = 0 $ and<br />
initial state $ u(t = 0) = 1$.<br />
<br />
An analytical solution for this domain is known, and we use it to evaluate or<br />
own solution. <br />
\begin{equation} <br />
u(\vec{p}, t) = \sum_{\substack{n=1 \\ n \text{<br />
odd}}}^\infty\sum_{\substack{m=1 \\ m \text{ odd}}}^\infty \frac{1}{\pi^2}<br />
\frac{16 a^2}{nm} \sin\left(\frac{\pi n}{a}p_x\right) \sin\left(\frac{\pi<br />
m}{a}p_y\right) e^{-\frac{\pi^2 (n^2+m^2)}{a^2}t} <br />
\end{equation}<br />
Because the solution is<br />
given in the series form, we only compare to the finite approximation, summing<br />
to $N = 100$ instead of infinity. Solution is on <xr id="fig:square_heat"/>.<br />
See the code for solving diffusion [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp here].<br />
<br />
<figure id="fig:square_heat"><br />
[[File:square_heat.png|thumb|upright=2|center|alt=A square of nodes coloured according to the solution(with smaller and larger node density)|<caption>A picture of our solution (with smaller and larger node density)</caption>]]<br />
</figure><br />
<br />
==Convergence with respect to number of nodes==<br />
<br />
<figure id="fig:node_convergence"><br />
[[File:node_convergence.png|thumb|upright=2|center|alt=Graph of errors with respect to number of nodes|<caption>Convergence with respect to number of nodes</caption>]]<br />
</figure><br />
<br />
We tested the method with a fixed time step of $ \Delta t = 1\cdot 10^{-5}$<br />
on a unit square ($a = 1$). Results are on <xr id="fig:node_convergence"/>. Monomial basis of $6$ monomials was used and $12$<br />
closest nodes counted as support for each node. After more than $250$ nodes of<br />
discretization in each dimension the method diverges, which is expected. The<br />
stability criterion for diffusion equation in two dimensions is $\Delta t \leq<br />
\frac{1}{4} \Delta x^2$, where $\Delta x$ is the spatial discretization<br />
step in one dimension. In our case, at $250$ nodes per side, the right hand side<br />
yields $\frac{1}{4}\cdot\frac{1}{250}\cdot\frac{1}{250} = 4\times 10^{-6}$,<br />
so our method is stable within the expected region.<br />
<br />
<figure id="fig:node_convergence_5"><br />
[[File:node_convergence_5.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for Gaussian and Monomial basis</caption>]]<br />
</figure><br />
<br />
On <xr id="fig:node_convergence_5"/> is another image of convergence, this time using monomial basis $\{1, x,<br />
y, x^2, y^2\}$ and Gaussian basis with discretization step $ \Delta t = 5<br />
\cdot 10^{-6}$ and 5 support nodes. Total node count was $N = 2500$.<br />
Error was calculated after $0.01$ time units have elapsed.<br />
<br />
<figure id="fig:node_convergence_comparison"><br />
[[File:explicit_implicit_diffusion_comparison.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for explicit and implicit mixed boundary condition (Neumann is enforced in two different ways).</caption>]]<br />
</figure><br />
<br />
Figure 3 shows convergence of 1/4 of the whole solution with Dirichlet boundary conditions (as above). The nodes were spanning the area $[0, 0.5]^2$. As expected, higher number of nodes yields lower max error up to some maximal number of nodes. The latter is lower for explicit than for implicit method.<br />
<br />
==Convergence with respect to number of time steps==<br />
<br />
<figure id="fig:timestep_convergence"><br />
[[File:timestep_convergence.png|thumb|upright=2|center|<caption>Convergence with respect to different timesteps</caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count with different time steps on the same<br />
domain as above. Results are on <xr id="fig:timestep_convergence"/>. For large time steps the method diverges, but once it starts<br />
converging the precision increases steadily as the time step decreases, until it<br />
hits its lower limit. This behaviour is expected. The error was calculated<br />
against the analytical solution above after $0.005$ units of time have passed. A<br />
monomial basis up to order $2$ inclusive ($m = 6$) was used and the support<br />
size was $n = 12$.<br />
<br />
==Using Gaussian basis==<br />
<br />
<figure id="fig:gauss_sigma_dependence"><br />
[[File:gauss_sigma_dependence.png|thumb|upright=2|center|<caption> Graph of error with respect to Gaussian basis $\sigma$ </caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count of $N = 2500$ with spatial step<br />
discretization of $\Delta x = \frac{1}{50}$. We used the Gaussian basis of<br />
$m = 5$ functions with support size $n = 13$. Error was calculated<br />
against analytical solution above after $0.01$ time units. A time step of<br />
$\Delta t = 10^{-5}$ was used. Results are on <xr id="fig:gauss_sigma_dependence"/><br />
<br />
As we can see from the graph, there exists an interval where the choice of<br />
$\sigma$ does not matter much, but outside of this interval, the method diverges<br />
rapidly. Care from the user side must be taken, to choose $\sigma$<br />
appropriately, with respect to plot above.<br />
<br />
== Convergence with respect to shape parameters ==<br />
<br />
Choice of shape of the basis function can have significant effect on solution of the system.<br />
We tested the space of feasible shape parameters on the diffusion equation with $\Delta x = \Delta y = 0.01$<br />
and $m = n = 9$. The weight function does not have much influence here, except on the condition <br />
number and any singular values that get cut off because of that.<br />
<br />
Results for Gaussian, MQ and IMQ basis functions are presented in <xr id="fig:shape_space_gau"/>, <xr id="fig:shape_space_mq"/> and <xr id="fig:shape_space_imq"/>.<br />
<br />
<figure id="fig:shape_space_gau"><br />
[[File:sigma_depedance_error_gau.png|thumb|left|800px|<caption>Error as a function of Gaussian basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_mq"><br />
[[File:sigma_depedance_error_imq.png|thumb|800px|<caption>Error as a function of MQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mq.png|thumb|upright=1|left|800px|<caption>Error as a function of IMQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mon.png|thumb|upright=1|800px|<caption>Error as a function of Gaussian weight's $\sigma_W$ for Monomials of order 2.</caption>]]<br />
</figure><br />
<br />
<div style="clear:both"></div><br />
<br />
==Solving in 3D==<br />
<br />
A 3-dimensional case on domain $[0, 1]^3$ was tested on $N = 20^3$<br />
nodes, making the discretization step $\Delta x = 0.05$. Support size of<br />
$n=10$ with $m=10$ Gaussian basis functions was used. Their<br />
normalization parameter was $\sigma = 60\Delta x$. A time step of $\Delta<br />
t = 10^{-5}$ and an explicit Euler method was used to calculate the solution<br />
of to $0.01$ time units. Resulting function is on <xr id="fig:diffusion3d"/>.<br />
<br />
<br />
<figure id="fig:diffusion3d"><br />
[[File:diffusion3d.png|thumb|upright =2|center|<caption>A $ 3 $-dimensional soulution with an explicit Euler method</caption>]]<br />
</figure><br />
<br />
==Solving Dirichlet with [http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classEngineMLS.html EngineMLS]==<br />
Example code using explicit stepping to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp example code]<br />
<br />
==Solving Dirichlet with MLSM operators==<br />
Example code using explicit stepping and MLSM operators to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion_mlsm_operators.cpp example code]<br />
<br />
==Solving mixed boundary conditions with MLSM operators==<br />
By using MLSM operators and utilizing its `[http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classmm_1_1MLSM.html MLSM::neumann]` method we can also solve<br />
partial diffusion equations. Solution is on <xr id="fig:quarter_diffusion"/>.<br />
<br />
<figure id="fig:quarter_diffusion"><br />
[[File:quarter_diffusion.png|400px|thumb|upright=2|center|<caption>Example of solving using Neumann boundary conditions </caption>]]<br />
</figure><br />
<br />
Example code showing the use of Neumann boundary conditions: [https://gitlab.com/e62Lab/e62numcodes/blob/master/test/mlsm_operators_test.cpp example code]<br />
<br />
We also support more -- interesting -- domains :) On <xr id="fig:mikimouse_heat"/> we see a domain in a shape of Miki mouse.<br />
<br />
<figure id="fig:mikimouse_heat"><br />
[[File:mikimouse_heat.png|400px|thumb|upright=2|center|<caption>Miki mouse domain </caption>]]<br />
</figure><br />
<br />
<!--<br />
==Convergence analysis for one dimensional diffusion equation==<br />
We now tried to solve the diffusion equation in one dimension <br />
$u_t = u_{xx}$ on the interval $[0,1]$<br />
with Dirichlet boundary conditions<br />
$u(0, t) = u(1,t) = 0 $<br />
and an initial state<br />
$u(x,0) = \sin(\pi x).$<br />
<br />
An analytical solution for this domain is given in closed form <br />
\begin{equation}<br />
u(x,t) = \sin(\pi x) \exp(-\pi^2 t),<br />
\end{equation}<br />
so we can use it to evaluate our approximation.<br />
<br />
We tested the method with a fixed time step of $\Delta t = 1\cdot 10^{-8}$ on a unit interval, that was discretized into $ N $ nodes. Monomial basis of $m$ monomials was used and $n$ closest nodes counted as support for each node. <br />
<br />
<strong> Note: </strong> the results are <strong>not confirmed </strong> and serve merely as an orientation. Errors may have occured. <br />
<br />
==Convergence with respect to number of nodes and support size==<br />
<br />
Convergence of the method was tested on different number of nodes $N = 15, 20, 25, \ldots, 200$. First we tested the method, when support size is equal to number of monomials in the basis ($n =m$), which we call local interpolation, for different values of $n = 3, 5, 7$. We then compared that to local approximation, where we take $n = 12$ supporting nodes and we only change number of monomials. We used gaussian weight function with $\sigma = 1\cdot \Delta x$. Error was calculated after $0.001$ time units have elapsed. We can see the results on <xr id="fig:diffusionConvergence1d"/>.<br />
<br />
<figure id="fig:diffusionConvergence1d"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional ddiffusion equation using monomial basis.</caption>]]<br />
</figure><br />
<br />
Larger is the support size, faster is convergence. This rule applies for local interpolation as well as approximation. But when number of nodes is too large, the method diverges faster for larger support sizes. Error for approximation is bigger than for interpolation (if we compare it for same basis size), but the method diverges with interpolation for smaller number of nodes.<br />
<br />
==Convergence with respect to weight function==<br />
<br />
We also checked the convergence for different weight functions. We tested it on number of nodes $N = 15, 20, 25, \ldots, 200$, with monomial basis with $m = 5$, support size $n = 12$ and gaussian weight function with $\sigma = w\cdot \Delta x$ for $w = 0.1$. Error was calculated after $0.001$ time units have elapsed with $\Delta t = 10^{-8}$. We can see the results on <xr id="fig:diffusionConvergence1dWeights"/>.<br />
<br />
<figure id="fig:diffusionConvergence1dWeights"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12_weights.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional diffusion equation using monomial basis with respect to weight function. </caption>]]<br />
</figure><br />
<br />
When the $\sigma$ is too small, other nodes on the support domain are not taken into consideration, so taking $n = 12$ or $n = 3$ would give the same results. However when $\sigma$ is too large, all the nodes on the support domain are contributing and there is no significant difference between different weight functions.<br />
--><br />
<br />
= Solving Electrostatics =<br />
<br />
This example is taken from the [http://www.freefem.org/ff++/ftp/freefem++doc.pdf FreeFem++ manual (page 235)].<br />
<br />
Assuming there is no current and the charge distribution is time independent, the electric field $\b{E}$ satisfies<br />
\begin{equation}\label{eq:electrostatics}<br />
\b{\nabla}\cdot\b{E} = \frac{\rho}{\epsilon}, \quad \b{\nabla} \times \b{E} = 0<br />
\end{equation}<br />
where $\rho$ is the charge density and $\epsilon$ is the permittivity. If we introduce an electrostatics potential $\phi$ such that<br />
\begin{equation}<br />
\b{E} = -\b{\nabla}\phi,<br />
\end{equation}<br />
we can insert it into the first equation in (\ref{eq:electrostatics}) resulting in Poisson's equation<br />
\begin{equation}\label{eq:poisson}<br />
\b{\nabla}^2 \phi = -\frac{\rho}{\epsilon}.<br />
\end{equation}<br />
In the absence of unpaired electric charge equation (\ref{eq:poisson}) becomes Laplace's equation<br />
\begin{equation}<br />
\b{\nabla}^2 \phi = 0<br />
\end{equation}<br />
<br />
We now solve this equation for a circular enclosure with two rectangular holes. The boundary of the circular enclosure is held at constant potential $0$ V. The two rectangular holes are held at constant potentials $+1$ V and $-1$ V, respectively. A coloured scatter plot is available in figure <xr id="fig:electro_statics"/>.<br />
<br />
<br />
<div class="toccolours mw-collapsible mw-collapsed"> '''Code for electrostatics problem'''<br />
<div class="mw-collapsible-content"><br />
<br />
<syntaxhighlight lang="c++" line><br />
// domain parameters, size, support and monomial order<br />
int domain_size = 3000;<br />
size_t n = 15; <br />
size_t m = 3; <br />
<br />
double radius = 5.0;<br />
double width = 0.3;<br />
double height = 3.0;<br />
<br />
// extra boundary labels<br />
int LEFT = -10;<br />
int RIGHT = -20;<br />
<br />
// build circular domain<br />
CircleDomain<Vec2d> domain({0,0},radius);<br />
domain.fillUniformInterior(domain_size);<br />
domain.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
<br />
// build left rectangle<br />
RectangleDomain<Vec2d> left({-2-width,-height},{-2+width,height});<br />
left.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
left.types[left.types == BOUNDARY] = LEFT;<br />
domain.subtract(left);<br />
<br />
// build right rectangle<br />
RectangleDomain<Vec2d> right({2-width,-height},{2+width,height});<br />
right.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
right.types[right.types == BOUNDARY] = RIGHT;<br />
domain.subtract(right);<br />
<br />
// relax domain and find supports<br />
domain.relax(50, 1e-2, 1.0, 3, 1000);<br />
domain.findSupport(n);<br />
<br />
// get interior and boundary ranges<br />
Range<int> interior = domain.types == INTERNAL;<br />
Range<int> boundary = domain.types == BOUNDARY;<br />
Range<int> leftrect = domain.types == LEFT;<br />
Range<int> rightrect = domain.types == RIGHT;<br />
<br />
// initialize unknown concentration field and right-hand side<br />
VecXd phi(domain.size(),1);<br />
VecXd RHS(domain.size(),1);<br />
<br />
// initialize interior values (this is an important step)<br />
RHS[interior] = 0.0;<br />
<br />
// set dirichlet boundary conditions<br />
RHS[boundary] = 0.0;<br />
RHS[leftrect] = -1.0;<br />
RHS[rightrect] = 1.0;<br />
<br />
// prepare shape functions and laplacian<br />
std::vector<Triplet<double>> l_coeff;<br />
for (auto& c : interior) {<br />
Range<Vec2d> supp_domain = domain.positions[domain.support[c]];<br />
EngineMLS<Vec2d, Monomials, Gaussians> MLS(m,supp_domain,pow(domain.characteristicDistance(),2));<br />
VecXd shape = MLS.getShapeAt(supp_domain[0], {{2, 0}}) +<br />
MLS.getShapeAt(supp_domain[0], {{0, 2}});<br />
for (size_t i = 0; i < supp_domain.size(); ++i) {<br />
l_coeff.emplace_back(c, domain.support[c][i], shape[i]);<br />
}<br />
}<br />
<br />
// prepare dirichlet boundaries<br />
std::vector<Triplet<double>> b_coeff;<br />
for (auto& c : domain.types < 0) {<br />
b_coeff.emplace_back(c, c, 1.0);<br />
}<br />
<br />
// prepare matrices<br />
SparseMatrix<double> laplace_m(domain.size(), domain.size());<br />
laplace_m.setFromTriplets(l_coeff.begin(), l_coeff.end());<br />
SparseMatrix<double> boundary_m(domain.size(), domain.size());<br />
boundary_m.setFromTriplets(b_coeff.begin(), b_coeff.end());<br />
<br />
// draw<br />
std::thread th([&] { draw2D(domain.positions, phi); });<br />
<br />
// initialize solver<br />
Eigen::BiCGSTAB<SparseMatrix<double>> solver;<br />
<br />
// join matrices together<br />
SparseMatrix<double> tmp = boundary_m + laplace_m;<br />
solver.compute(tmp);<br />
<br />
// solve system of equations<br />
phi = solver.solve(RHS);<br />
std::cout << "#iterations: " << solver.iterations() << std::endl;<br />
std::cout << "estimated error: " << solver.error() << std::endl;<br />
<br />
// end drawing<br />
th.join(); <br />
</syntaxhighlight><br />
<br />
</div><br />
</div><br />
<br />
<br />
<figure id="fig:electro_statics"><br />
[[File:electrostatics.png|thumb|upright=2|center|<caption> Simple electrostatics problem solved by implicit method with 15 node monomial supports and gaussian weight functions. </caption>]]<br />
</figure><br />
<br />
= Convection-diffusion =<br />
<br />
The following example problem is adapted from [http://servidor.demec.ufpr.br/CFD/bibliografia/Ferziger%20Peric%20-%20Computational%20Methods%20for%20Fluid%20Dynamics,%203rd%20Ed%20-%202002.pdf Ferziger & Perič (2002) (pages 82-86)].<br />
<br />
<figure id="fig:fvm_problem"><br />
[[File:Screenshot_2016-11-28_11-05-03.png|thumb|400px|upright=2|centre|<caption> Geometry and boundary conditions for the scalar transport in a stagnation point flow. </caption>]]<br />
</figure><br />
<br />
We now consider the problem of a scalar quantity transported in a known velocity field $\b{u}$. The velocity field is given by<br />
\[\b{u} = (u_x,u_y) = (x, -y),\]<br />
which represents the flow near a stagnation point $(x,y) = (0,0)$. The streamlines $xy=\mathrm{const.}$ and change direction with respect to the Cartesian grid. The convection-diffusion equation (also known as the scalar transport equation) to be solved is<br />
\[\alpha \b{\nabla}^2 \phi - \b{u}\cdot \b{\nabla}\phi = 0\]<br />
with the following boundary conditions:<br />
* $\phi = 0$ along the north (inlet) boundary;<br />
* linear variation $\phi = (1 - y)$ along the west boundary ($x = 0$);<br />
* symmetry condition on the south boundary;<br />
* zero gradient at the east (outlet) boundary.<br />
The parameter $\alpha$ is the diffusivity for the scalar quantity $\phi$. The problem and boundary conditions are schematically represented in the <xr id="fig:fvm_problem"/>. The solutions obtained by Ferziger and Perić are shown in <xr id="fig:fvm_problem_solution"/>.<br />
<br />
<figure id="fig:fvm_problem_solution"><br />
[[File:Screenshot_2016-11-28_11-07-44.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right). </caption>]]<br />
</figure><br />
<br />
The meshless solution is obtained on regular point arrangements of sizes 100 × 100 and 300 × 300, for the 2 cases of diffusivity $\alpha = 0.01$ and $\alpha = 0.001$, respectively. For comparison we have plotted the contour lines at the same intervals as the original authors of this example (see <xr id="fig:fvm_problem_meshless"/>). The basis functions are the 5 monomials $1$, $x$, $y$, $x^2$, $y^2$.<br />
<br />
<figure id="fig:fvm_problem_meshless"><br />
[[File:fvm_convection.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right), obtaines with meshless solver. </caption>]]<br />
</figure></div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=Analysis_of_MLSM_performance&diff=1125Analysis of MLSM performance2017-04-24T09:00:57Z<p>JLapajne: /* Convergence with respect to number of nodes */</p>
<hr />
<div>=Solving Diffusion equation=<br />
For starters, we can solve simple Diffusion equation<br />
$ \nabla^2 u = \frac{\partial u}{\partial t} $.<br />
<br />
We solved the equation on a square $\Omega = [0, a] \times [0, a]$ with<br />
Dirichlet boundary conditions $ \left. u\right|_{\partial \Omega} = 0 $ and<br />
initial state $ u(t = 0) = 1$.<br />
<br />
An analytical solution for this domain is known, and we use it to evaluate or<br />
own solution. <br />
\begin{equation} <br />
u(\vec{p}, t) = \sum_{\substack{n=1 \\ n \text{<br />
odd}}}^\infty\sum_{\substack{m=1 \\ m \text{ odd}}}^\infty \frac{1}{\pi^2}<br />
\frac{16 a^2}{nm} \sin\left(\frac{\pi n}{a}p_x\right) \sin\left(\frac{\pi<br />
m}{a}p_y\right) e^{-\frac{\pi^2 (n^2+m^2)}{a^2}t} <br />
\end{equation}<br />
Because the solution is<br />
given in the series form, we only compare to the finite approximation, summing<br />
to $N = 100$ instead of infinity. Solution is on <xr id="fig:square_heat"/>.<br />
See the code for solving diffusion [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp here].<br />
<br />
<figure id="fig:square_heat"><br />
[[File:square_heat.png|thumb|upright=2|center|alt=A square of nodes coloured according to the solution(with smaller and larger node density)|<caption>A picture of our solution (with smaller and larger node density)</caption>]]<br />
</figure><br />
<br />
==Convergence with respect to number of nodes==<br />
<br />
<figure id="fig:node_convergence"><br />
[[File:node_convergence.png|thumb|upright=2|center|alt=Graph of errors with respect to number of nodes|<caption>Convergence with respect to number of nodes</caption>]]<br />
</figure><br />
<br />
We tested the method with a fixed time step of $ \Delta t = 1\cdot 10^{-5}$<br />
on a unit square ($a = 1$). Results are on <xr id="fig:node_convergence"/>. Monomial basis of $6$ monomials was used and $12$<br />
closest nodes counted as support for each node. After more than $250$ nodes of<br />
discretization in each dimension the method diverges, which is expected. The<br />
stability criterion for diffusion equation in two dimensions is $\Delta t \leq<br />
\frac{1}{4} \Delta x^2$, where $\Delta x$ is the spatial discretization<br />
step in one dimension. In our case, at $250$ nodes per side, the right hand side<br />
yields $\frac{1}{4}\cdot\frac{1}{250}\cdot\frac{1}{250} = 4\times 10^{-6}$,<br />
so our method is stable within the expected region.<br />
<br />
<figure id="fig:node_convergence_5"><br />
[[File:node_convergence_5.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for Gaussian and Monomial basis</caption>]]<br />
</figure><br />
<br />
On <xr id="fig:node_convergence_5"/> is another image of convergence, this time using monomial basis $\{1, x,<br />
y, x^2, y^2\}$ and Gaussian basis with discretization step $ \Delta t = 5<br />
\cdot 10^{-6}$ and 5 support nodes. Total node count was $N = 2500$.<br />
Error was calculated after $0.01$ time units have elapsed.<br />
<br />
<figure id="fig:node_convergence_5"><br />
[[File:explicit_implicit_diffusion_comparison.png|thumb|upright=2|center|<caption>Convergence with respect to number of nodes for explicit and implicit mixed boundary condition (Neumann is enforced in two different ways).</caption>]]<br />
</figure><br />
<br />
Figure 3 shows convergence of 1/4 of the whole solution with Dirichlet boundary conditions (as above). The nodes were spanning the area $[0, 0.5]^2$. As expected, higher number of nodes yields lower max error up to some maximal number of nodes. The latter is lower for explicit than for implicit method.<br />
<br />
==Convergence with respect to number of time steps==<br />
<br />
<figure id="fig:timestep_convergence"><br />
[[File:timestep_convergence.png|thumb|upright=2|center|<caption>Convergence with respect to different timesteps</caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count with different time steps on the same<br />
domain as above. Results are on <xr id="fig:timestep_convergence"/>. For large time steps the method diverges, but once it starts<br />
converging the precision increases steadily as the time step decreases, until it<br />
hits its lower limit. This behaviour is expected. The error was calculated<br />
against the analytical solution above after $0.005$ units of time have passed. A<br />
monomial basis up to order $2$ inclusive ($m = 6$) was used and the support<br />
size was $n = 12$.<br />
<br />
==Using Gaussian basis==<br />
<br />
<figure id="fig:gauss_sigma_dependence"><br />
[[File:gauss_sigma_dependence.png|thumb|upright=2|center|<caption> Graph of error with respect to Gaussian basis $\sigma$ </caption>]]<br />
</figure><br />
<br />
We tested the method on a fixed node count of $N = 2500$ with spatial step<br />
discretization of $\Delta x = \frac{1}{50}$. We used the Gaussian basis of<br />
$m = 5$ functions with support size $n = 13$. Error was calculated<br />
against analytical solution above after $0.01$ time units. A time step of<br />
$\Delta t = 10^{-5}$ was used. Results are on <xr id="fig:gauss_sigma_dependence"/><br />
<br />
As we can see from the graph, there exists an interval where the choice of<br />
$\sigma$ does not matter much, but outside of this interval, the method diverges<br />
rapidly. Care from the user side must be taken, to choose $\sigma$<br />
appropriately, with respect to plot above.<br />
<br />
== Convergence with respect to shape parameters ==<br />
<br />
Choice of shape of the basis function can have significant effect on solution of the system.<br />
We tested the space of feasible shape parameters on the diffusion equation with $\Delta x = \Delta y = 0.01$<br />
and $m = n = 9$. The weight function does not have much influence here, except on the condition <br />
number and any singular values that get cut off because of that.<br />
<br />
Results for Gaussian, MQ and IMQ basis functions are presented in <xr id="fig:shape_space_gau"/>, <xr id="fig:shape_space_mq"/> and <xr id="fig:shape_space_imq"/>.<br />
<br />
<figure id="fig:shape_space_gau"><br />
[[File:sigma_depedance_error_gau.png|thumb|left|800px|<caption>Error as a function of Gaussian basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_mq"><br />
[[File:sigma_depedance_error_imq.png|thumb|800px|<caption>Error as a function of MQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mq.png|thumb|upright=1|left|800px|<caption>Error as a function of IMQ basis' $\sigma_B$ and Gaussian weight's $\sigma_W$</caption>]]<br />
</figure><br />
<figure id="fig:shape_space_imq"><br />
[[File:sigma_depedance_error_mon.png|thumb|upright=1|800px|<caption>Error as a function of Gaussian weight's $\sigma_W$ for Monomials of order 2.</caption>]]<br />
</figure><br />
<br />
<div style="clear:both"></div><br />
<br />
==Solving in 3D==<br />
<br />
A 3-dimensional case on domain $[0, 1]^3$ was tested on $N = 20^3$<br />
nodes, making the discretization step $\Delta x = 0.05$. Support size of<br />
$n=10$ with $m=10$ Gaussian basis functions was used. Their<br />
normalization parameter was $\sigma = 60\Delta x$. A time step of $\Delta<br />
t = 10^{-5}$ and an explicit Euler method was used to calculate the solution<br />
of to $0.01$ time units. Resulting function is on <xr id="fig:diffusion3d"/>.<br />
<br />
<br />
<figure id="fig:diffusion3d"><br />
[[File:diffusion3d.png|thumb|upright =2|center|<caption>A $ 3 $-dimensional soulution with an explicit Euler method</caption>]]<br />
</figure><br />
<br />
==Solving Dirichlet with [http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classEngineMLS.html EngineMLS]==<br />
Example code using explicit stepping to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion.cpp example code]<br />
<br />
==Solving Dirichlet with MLSM operators==<br />
Example code using explicit stepping and MLSM operators to reproduce the thermal images from the beginning: [https://gitlab.com/e62Lab/e62numcodes/blob/master/examples/diffusion/diffusion_mlsm_operators.cpp example code]<br />
<br />
==Solving mixed boundary conditions with MLSM operators==<br />
By using MLSM operators and utilizing its `[http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/classmm_1_1MLSM.html MLSM::neumann]` method we can also solve<br />
partial diffusion equations. Solution is on <xr id="fig:quarter_diffusion"/>.<br />
<br />
<figure id="fig:quarter_diffusion"><br />
[[File:quarter_diffusion.png|400px|thumb|upright=2|center|<caption>Example of solving using Neumann boundary conditions </caption>]]<br />
</figure><br />
<br />
Example code showing the use of Neumann boundary conditions: [https://gitlab.com/e62Lab/e62numcodes/blob/master/test/mlsm_operators_test.cpp example code]<br />
<br />
We also support more -- interesting -- domains :) On <xr id="fig:mikimouse_heat"/> we see a domain in a shape of Miki mouse.<br />
<br />
<figure id="fig:mikimouse_heat"><br />
[[File:mikimouse_heat.png|400px|thumb|upright=2|center|<caption>Miki mouse domain </caption>]]<br />
</figure><br />
<br />
<!--<br />
==Convergence analysis for one dimensional diffusion equation==<br />
We now tried to solve the diffusion equation in one dimension <br />
$u_t = u_{xx}$ on the interval $[0,1]$<br />
with Dirichlet boundary conditions<br />
$u(0, t) = u(1,t) = 0 $<br />
and an initial state<br />
$u(x,0) = \sin(\pi x).$<br />
<br />
An analytical solution for this domain is given in closed form <br />
\begin{equation}<br />
u(x,t) = \sin(\pi x) \exp(-\pi^2 t),<br />
\end{equation}<br />
so we can use it to evaluate our approximation.<br />
<br />
We tested the method with a fixed time step of $\Delta t = 1\cdot 10^{-8}$ on a unit interval, that was discretized into $ N $ nodes. Monomial basis of $m$ monomials was used and $n$ closest nodes counted as support for each node. <br />
<br />
<strong> Note: </strong> the results are <strong>not confirmed </strong> and serve merely as an orientation. Errors may have occured. <br />
<br />
==Convergence with respect to number of nodes and support size==<br />
<br />
Convergence of the method was tested on different number of nodes $N = 15, 20, 25, \ldots, 200$. First we tested the method, when support size is equal to number of monomials in the basis ($n =m$), which we call local interpolation, for different values of $n = 3, 5, 7$. We then compared that to local approximation, where we take $n = 12$ supporting nodes and we only change number of monomials. We used gaussian weight function with $\sigma = 1\cdot \Delta x$. Error was calculated after $0.001$ time units have elapsed. We can see the results on <xr id="fig:diffusionConvergence1d"/>.<br />
<br />
<figure id="fig:diffusionConvergence1d"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional ddiffusion equation using monomial basis.</caption>]]<br />
</figure><br />
<br />
Larger is the support size, faster is convergence. This rule applies for local interpolation as well as approximation. But when number of nodes is too large, the method diverges faster for larger support sizes. Error for approximation is bigger than for interpolation (if we compare it for same basis size), but the method diverges with interpolation for smaller number of nodes.<br />
<br />
==Convergence with respect to weight function==<br />
<br />
We also checked the convergence for different weight functions. We tested it on number of nodes $N = 15, 20, 25, \ldots, 200$, with monomial basis with $m = 5$, support size $n = 12$ and gaussian weight function with $\sigma = w\cdot \Delta x$ for $w = 0.1$. Error was calculated after $0.001$ time units have elapsed with $\Delta t = 10^{-8}$. We can see the results on <xr id="fig:diffusionConvergence1dWeights"/>.<br />
<br />
<figure id="fig:diffusionConvergence1dWeights"><br />
[[File:monomials_convergnece_dt_e-8_t_0-001_approxSupp_12_weights.png|thumb|upright=2|center|<caption> Convergence analysis for $ 1 $-dimensional diffusion equation using monomial basis with respect to weight function. </caption>]]<br />
</figure><br />
<br />
When the $\sigma$ is too small, other nodes on the support domain are not taken into consideration, so taking $n = 12$ or $n = 3$ would give the same results. However when $\sigma$ is too large, all the nodes on the support domain are contributing and there is no significant difference between different weight functions.<br />
--><br />
<br />
= Solving Electrostatics =<br />
<br />
This example is taken from the [http://www.freefem.org/ff++/ftp/freefem++doc.pdf FreeFem++ manual (page 235)].<br />
<br />
Assuming there is no current and the charge distribution is time independent, the electric field $\b{E}$ satisfies<br />
\begin{equation}\label{eq:electrostatics}<br />
\b{\nabla}\cdot\b{E} = \frac{\rho}{\epsilon}, \quad \b{\nabla} \times \b{E} = 0<br />
\end{equation}<br />
where $\rho$ is the charge density and $\epsilon$ is the permittivity. If we introduce an electrostatics potential $\phi$ such that<br />
\begin{equation}<br />
\b{E} = -\b{\nabla}\phi,<br />
\end{equation}<br />
we can insert it into the first equation in (\ref{eq:electrostatics}) resulting in Poisson's equation<br />
\begin{equation}\label{eq:poisson}<br />
\b{\nabla}^2 \phi = -\frac{\rho}{\epsilon}.<br />
\end{equation}<br />
In the absence of unpaired electric charge equation (\ref{eq:poisson}) becomes Laplace's equation<br />
\begin{equation}<br />
\b{\nabla}^2 \phi = 0<br />
\end{equation}<br />
<br />
We now solve this equation for a circular enclosure with two rectangular holes. The boundary of the circular enclosure is held at constant potential $0$ V. The two rectangular holes are held at constant potentials $+1$ V and $-1$ V, respectively. A coloured scatter plot is available in figure <xr id="fig:electro_statics"/>.<br />
<br />
<br />
<div class="toccolours mw-collapsible mw-collapsed"> '''Code for electrostatics problem'''<br />
<div class="mw-collapsible-content"><br />
<br />
<syntaxhighlight lang="c++" line><br />
// domain parameters, size, support and monomial order<br />
int domain_size = 3000;<br />
size_t n = 15; <br />
size_t m = 3; <br />
<br />
double radius = 5.0;<br />
double width = 0.3;<br />
double height = 3.0;<br />
<br />
// extra boundary labels<br />
int LEFT = -10;<br />
int RIGHT = -20;<br />
<br />
// build circular domain<br />
CircleDomain<Vec2d> domain({0,0},radius);<br />
domain.fillUniformInterior(domain_size);<br />
domain.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
<br />
// build left rectangle<br />
RectangleDomain<Vec2d> left({-2-width,-height},{-2+width,height});<br />
left.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
left.types[left.types == BOUNDARY] = LEFT;<br />
domain.subtract(left);<br />
<br />
// build right rectangle<br />
RectangleDomain<Vec2d> right({2-width,-height},{2+width,height});<br />
right.fillUniformBoundaryWithStep(domain.characteristicDistance());<br />
right.types[right.types == BOUNDARY] = RIGHT;<br />
domain.subtract(right);<br />
<br />
// relax domain and find supports<br />
domain.relax(50, 1e-2, 1.0, 3, 1000);<br />
domain.findSupport(n);<br />
<br />
// get interior and boundary ranges<br />
Range<int> interior = domain.types == INTERNAL;<br />
Range<int> boundary = domain.types == BOUNDARY;<br />
Range<int> leftrect = domain.types == LEFT;<br />
Range<int> rightrect = domain.types == RIGHT;<br />
<br />
// initialize unknown concentration field and right-hand side<br />
VecXd phi(domain.size(),1);<br />
VecXd RHS(domain.size(),1);<br />
<br />
// initialize interior values (this is an important step)<br />
RHS[interior] = 0.0;<br />
<br />
// set dirichlet boundary conditions<br />
RHS[boundary] = 0.0;<br />
RHS[leftrect] = -1.0;<br />
RHS[rightrect] = 1.0;<br />
<br />
// prepare shape functions and laplacian<br />
std::vector<Triplet<double>> l_coeff;<br />
for (auto& c : interior) {<br />
Range<Vec2d> supp_domain = domain.positions[domain.support[c]];<br />
EngineMLS<Vec2d, Monomials, Gaussians> MLS(m,supp_domain,pow(domain.characteristicDistance(),2));<br />
VecXd shape = MLS.getShapeAt(supp_domain[0], {{2, 0}}) +<br />
MLS.getShapeAt(supp_domain[0], {{0, 2}});<br />
for (size_t i = 0; i < supp_domain.size(); ++i) {<br />
l_coeff.emplace_back(c, domain.support[c][i], shape[i]);<br />
}<br />
}<br />
<br />
// prepare dirichlet boundaries<br />
std::vector<Triplet<double>> b_coeff;<br />
for (auto& c : domain.types < 0) {<br />
b_coeff.emplace_back(c, c, 1.0);<br />
}<br />
<br />
// prepare matrices<br />
SparseMatrix<double> laplace_m(domain.size(), domain.size());<br />
laplace_m.setFromTriplets(l_coeff.begin(), l_coeff.end());<br />
SparseMatrix<double> boundary_m(domain.size(), domain.size());<br />
boundary_m.setFromTriplets(b_coeff.begin(), b_coeff.end());<br />
<br />
// draw<br />
std::thread th([&] { draw2D(domain.positions, phi); });<br />
<br />
// initialize solver<br />
Eigen::BiCGSTAB<SparseMatrix<double>> solver;<br />
<br />
// join matrices together<br />
SparseMatrix<double> tmp = boundary_m + laplace_m;<br />
solver.compute(tmp);<br />
<br />
// solve system of equations<br />
phi = solver.solve(RHS);<br />
std::cout << "#iterations: " << solver.iterations() << std::endl;<br />
std::cout << "estimated error: " << solver.error() << std::endl;<br />
<br />
// end drawing<br />
th.join(); <br />
</syntaxhighlight><br />
<br />
</div><br />
</div><br />
<br />
<br />
<figure id="fig:electro_statics"><br />
[[File:electrostatics.png|thumb|upright=2|center|<caption> Simple electrostatics problem solved by implicit method with 15 node monomial supports and gaussian weight functions. </caption>]]<br />
</figure><br />
<br />
= Convection-diffusion =<br />
<br />
The following example problem is adapted from [http://servidor.demec.ufpr.br/CFD/bibliografia/Ferziger%20Peric%20-%20Computational%20Methods%20for%20Fluid%20Dynamics,%203rd%20Ed%20-%202002.pdf Ferziger & Perič (2002) (pages 82-86)].<br />
<br />
<figure id="fig:fvm_problem"><br />
[[File:Screenshot_2016-11-28_11-05-03.png|thumb|400px|upright=2|centre|<caption> Geometry and boundary conditions for the scalar transport in a stagnation point flow. </caption>]]<br />
</figure><br />
<br />
We now consider the problem of a scalar quantity transported in a known velocity field $\b{u}$. The velocity field is given by<br />
\[\b{u} = (u_x,u_y) = (x, -y),\]<br />
which represents the flow near a stagnation point $(x,y) = (0,0)$. The streamlines $xy=\mathrm{const.}$ and change direction with respect to the Cartesian grid. The convection-diffusion equation (also known as the scalar transport equation) to be solved is<br />
\[\alpha \b{\nabla}^2 \phi - \b{u}\cdot \b{\nabla}\phi = 0\]<br />
with the following boundary conditions:<br />
* $\phi = 0$ along the north (inlet) boundary;<br />
* linear variation $\phi = (1 - y)$ along the west boundary ($x = 0$);<br />
* symmetry condition on the south boundary;<br />
* zero gradient at the east (outlet) boundary.<br />
The parameter $\alpha$ is the diffusivity for the scalar quantity $\phi$. The problem and boundary conditions are schematically represented in the <xr id="fig:fvm_problem"/>. The solutions obtained by Ferziger and Perić are shown in <xr id="fig:fvm_problem_solution"/>.<br />
<br />
<figure id="fig:fvm_problem_solution"><br />
[[File:Screenshot_2016-11-28_11-07-44.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right). </caption>]]<br />
</figure><br />
<br />
The meshless solution is obtained on regular point arrangements of sizes 100 × 100 and 300 × 300, for the 2 cases of diffusivity $\alpha = 0.01$ and $\alpha = 0.001$, respectively. For comparison we have plotted the contour lines at the same intervals as the original authors of this example (see <xr id="fig:fvm_problem_meshless"/>). The basis functions are the 5 monomials $1$, $x$, $y$, $x^2$, $y^2$.<br />
<br />
<figure id="fig:fvm_problem_meshless"><br />
[[File:fvm_convection.png|thumb|600px|center|upright=2|<caption> Isolines of $\phi$, from 0.05 to 0.95 with step 0.1 (top to bottom), for $\alpha = 0.01$ (left) and $\alpha = 0.001$ (right), obtaines with meshless solver. </caption>]]<br />
</figure></div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1013K-d tree2017-03-17T08:24:40Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster.<br />
<br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
Nearest neighbor search is a $logn$ procedure, although in practice, especially in higher dimensions, it's closer to o(n). There is a nice description of the algorithm on [https://en.wikipedia.org/wiki/K-d_tree wikipedia] website. We tested ANN and Libssrckdtree implementations. In addition we tested also cover tree nearest neighbor search.<br />
<br />
[[File:nearest_neighbor_search_time.jpg|800px]]<br />
<br />
Note that cover tree implementation is not the best available, but it is the only one to have insert/remove and is written in c++.<br />
<br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.<br />
*k-d tree. (n.d.). In Wikipedia. Retrieved March 17, 2017, from https://en.wikipedia.org/wiki/K-d_tree</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1012K-d tree2017-03-17T08:23:03Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster.<br />
<br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
Nearest neighbor search is a $logn$ procedure, although in practice, especially in higher dimensions, it's closer to o(n). There is a nice description of the algorithm on [https://en.wikipedia.org/wiki/K-d_tree wikipedia] website. We tested ANN and Libssrckdtree implementations. In addition we tested also cover tree nearest neighbor search.<br />
<br />
[[File:nearest_neighbor_search_time.jpg|800px]]<br />
<br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.<br />
*k-d tree. (n.d.). In Wikipedia. Retrieved March 17, 2017, from https://en.wikipedia.org/wiki/K-d_tree</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1011K-d tree2017-03-17T08:20:22Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster.<br />
<br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
Nearest neighbor search is a $logn$ procedure, although in practice, especially in higher dimensions, it's closer to o(n). There is a nice description of the algorithm on [https://en.wikipedia.org/wiki/K-d_tree wikipedia] website. We tested ANN and Libssrckdtree implementations. In addition we tested also cover tree nearest neighbor search.<br />
<br />
[[File:nearest_neighbor_search_time.jpg|800px]]<br />
<br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1010K-d tree2017-03-17T08:19:26Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster.<br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
Nearest neighbor search is a $logn$ procedure, although in practice, especially in higher dimensions, it's closer to o(n). There is a nice description of the algorithm on [https://en.wikipedia.org/wiki/K-d_tree wikipedia] website. We tested ANN and Libssrckdtree implementations. In addition we tested also cover tree nearest neighbor search.<br />
<br />
[[File:nearest_neighbor_search_time.jpg|800px]]<br />
<br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1009K-d tree2017-03-17T08:18:20Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
Nearest neighbor search is a $logn$ procedure, although in practice, especially in higher dimensions, it's closer to o(n). There is a nice description of the algorithm on [https://en.wikipedia.org/wiki/K-d_tree wikipedia] website. We tested ANN and Libssrckdtree implementations. In addition we tested also cover tree nearest neighbor search.<br />
<br />
[[File:nearest_neighbor_search_time.jpg|800px]]<br />
<br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1007K-d tree2017-03-17T08:13:46Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
Nearest neighbor search is a $logn$ procedure. There is a nice description of the algorithm on [https://en.wikipedia.org/wiki/K-d_tree wikipedia] website <br />
<br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1006K-d tree2017-03-17T08:01:28Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:Tree_building_time.jpg]]<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1004K-d tree2017-03-17T07:58:45Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:cas_gradnje.jpg]]<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1003K-d tree2017-03-17T07:57:47Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:cas_gradnje.pdf]]<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1002K-d tree2017-03-17T07:56:25Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We also included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
[[File:/home/jure/delo_ijs_e6/libssrckdtree_test/nearest_neighbors_cas\cas_gradnje.pdf������������������]]<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1001K-d tree2017-03-17T07:53:35Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding (consult literature). Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries. We included a cover tree library, which can be found [https://github.com/DNCrane/Cover-Tree here].<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=1000K-d tree2017-03-17T07:49:25Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding. Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
We tested build times for Ann and Libssrckdtree libraries.<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=999K-d tree2017-03-17T07:46:46Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding. Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o($log n$). Inserting n nodes thus again results in $n log n$ time complexity.<br />
<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=998K-d tree2017-03-17T07:45:45Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The procedure described above is used to construct a tree from a collection of given points and is used by ANN library. The tree can also be constructed by inserting points one by one. Such procedure is used by Libssrckdtree library. Time complexity of both procedures is $o(nlogn)$. The most time consuming part of the former algorithm is median point selection. It's time complexity is o(n). This can be improved by imploying some smarter schemes than median finding. Time complexity of the second algorithm can quickly be calculated from the average time need to insert a node into the tree - it scales o(log n). Inserting n nodes thus again results in $nlogn$ time complexity.<br />
<br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=997K-d tree2017-03-17T07:18:03Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the node in a tree. At depth $k$, the k-th axis (dimenison) is used to divide the remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level, we look for median point by looking at x coordiantes of all points. We put the median point into the root leaf. Remaining points are put into the left and right subtrees. We put the points which have smaller x coorinate into the left subtree and points which have larger x coorinate into the right subtree (or vice versa). We construct the two children nodes of the root from the two subgroups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. We continue this way until we run out of points. The axes used for median search are periodically repeating.<br />
<br />
The above procedure is the most commond <br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=950K-d tree2017-03-15T11:52:21Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the tree. At depth $k$, the k-th axis (dimenison) is used to divide remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level we look for median point by looking at x coordiantes of points. We put the median points into the root leaf. Remaining points are put into the left and right subtrees. In the left subtree we put the points which have smaller x coorinate and in the right subtree the points which have larger x coordinate (or vice versa). We construct the two children nodes of the root from the two groups we just created. The procedure is the same, except that we divide the space according to y or z coordinate. <br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajnehttp://e6.ijs.si/medusa/wiki/index.php?title=K-d_tree&diff=949K-d tree2017-03-15T11:47:56Z<p>JLapajne: </p>
<hr />
<div>In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space. <br />
In unstructured mesh based approaches, support domains can be determined from the mesh topology. However, in meshless methods, the nearest neighboring nodes in $\Omega_S$ are determined with various algorithms and specialized data structures, we use kD Tree.<br />
<br />
Kd tree is a binary tree which uses hyperplane space partitioning to divide points into subgroups (subtrees). The plane is chosen based on the depth of the tree. At depth $k$, the k-th axis (dimenison) is used to divide remaining points. Example: in 3 dimensions we have x,y and z axes. At root level we decide to divide points according to their $x$ values (we could also do it according to y or z values). Usual criterion for divison of points is median. Thus, at root level we look for median point by looking at x coordiantes of points. We put the median points into the root leaf. Remaining points are put into the left and right subtrees. We put in the left subtree the points which have smaller x coorinate and in the right subtree the points which have larger x coordinate (or vice versa). <br />
<br />
The input to Algorithm is a list of $(n_S +1)$ nearest nodes to ${\bf x}$. A naive implementation of finding these nodes could become quite complex if implemented with classical sorting of all nodal distances. Regarding the number of nodes $N$, it approaches the quadratic computational complexity. However, by using efficient data structures, e.g. quadtree, R tree or kD tree (2D tree for 2D domains), the problem becomes tractable. <br />
The strategy is to build the data structure only once, before the solution procedure. During the solution, the support nodes of the desired points will then be found much faster. <br />
<br />
<!--<br />
Let us illustrate the whole procedure with a simple example of a 2D tree for eleven nodes, with node numbers and coordinates listed in the first and the second column of Table \ref{tab:kD_nodes} and attached to the corresponding dots on the unit square in the left part of Figure \ref{fig.2D_treegraph}. <br />
<br />
{\begin{table}[h]<br />
\centering<br />
\begin{tabular}{|c|c|c|c|c|c|}<br />
\hline<br />
Node number&Unsorted & Sorted by $x$ & Sorted by $y$ & Sorted by $x$ & Bucket\\ \hline<br />
1&(0,0)&(0,0)&(0,0)&(0,0)& \bf (0,0) \\ \hline<br />
2&(0.6,0)&(0,0.4)&(0,0.4)&\bf (0,0.4)& \\ \hline<br />
3&(1,0)&(0,1)&\bf (0.24,0.6)& & \\ \hline<br />
4&(0,0.4)&(0.24,0.6)&(0,1)&(0,1)& \bf (0,1) \\ \hline<br />
5&(0.6,0.3)&(0.47,1)&(0.47,1)&\bf (0.47,1)&\\ \hline<br />
6&(1,0.5)&\bf (0.6,0)& & &\\ \hline<br />
7&(0.24,0.6)&(0.6,0.3)&(1,0)&(0.6,0.3)&\bf (0.6,0.3)\\ \hline<br />
8&(0.76,0.8)&(0.76,0.8)&(0.6,0.3)&\bf (1,0)&\\ \hline<br />
9&(0,1)&(1,0)&\bf (1,0.5)& &\\ \hline<br />
10&(0.47,1)&(1,0.5)&(0.76,0.8)&(0.76,0.8)&\bf (0.76,0.8)\\ \hline<br />
11&(1,1)&(1,1)&(1,1)&\bf (1,1)&\\ \hline<br />
\end{tabular}<br />
\caption{The list of eleven nodes (1st column) determined with coordinates (2nd column), after sorting by $x$ (3rd column), after sorting sub lists by $y$ (4th column) and after sorting sub sub lists by $x$ again (5th column). The nodes nearest to medians are shown in bold.}<br />
\label{tab:kD_nodes}<br />
\end{table}<br />
<br />
In the first step of 2D tree construction, the list of nodes is sorted by their $x$ coordinate, which is shown in third column of Table \ref{tab:kD_nodes}. Then a node with median coordinate, $x = 0.6$ in our case (shown in bold), is selected as the root of the first level of the 2D tree. If there is more than one such node, any one can be selected. The sorted set in column 3 is split into two parts, the one for $x$ below the median, i.e. $x < 0.6$, and the one for $x$ above or equal the median, i.e. $x\geq 0.6$. The two sub sets of nodes are shown in the left part of Figure \ref{fig.2D_treegraph} within two distinct rectangles, and on the right side of Figure \ref{fig.2D_treegraph} as the left and the right part of the 2D tree. <br />
<br />
In the second step, the two sub lists of nodes are sorted by their $y$ coordinate, which is shown in fourth column of Table \ref{tab:kD_nodes}. The median coordinates $y$ are $0.6$ and $0.5$, respectively. The corresponding nodes $(0.24,0.6)$ and $(1,0.5)$ are taken as roots for the second level of 2D tree and are shown in bold and used to further split the tree. The resulting four sub sub sets of nodes are shown in the right side of Figure \ref{fig.2D_treegraph} as nodes on the lower two levels of the 2D tree. <br />
<br />
Finally, the sub sub lists are sorted again by their $x$ coordinate, with the result shown in fifth column. Four roots are obtained with the coordinate $x$ nearest to medians, namely the nodes $(0,0.4)$, $(0.47,1)$, $(1,0)$, and $(1,1)$. The remaining nodes of the last level of the 2D tree, also termed the bucket, are<br />
$(0,0)$, $(0,1)$, $(0.6,0.3)$, and $(0.76,0.8)$.<br />
In practical cases, the refinement of the tree stops sooner, when its leaves are represented by list of several nodes, because such a fine grained distribution of leaves as in the presented example is often not beneficial from the computational efficiency point of view.<br />
--><br />
[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]<br />
<br />
== References ==<br />
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321;<br />
* Trobec R., Kosec G., Parallel scientific computing : theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.</div>JLapajne