Medusa: Coordinate Free Mehless Method implementation - User contributions [en]

Scattering from an infinite cylinder

2026-01-22T16:24:37Z

AndrejKP: Fix the boundary condition arising from H fields (no epsilon there) and the incident wave form. Add solution reference.

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad D
\label{eq:outer}
\end{align}
The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =E_z^{inc} \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial D.
\end{equation}
The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{-ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

<syntaxhighlight lang="c++" line>
// Construct system
double l_in = epsr*std::pow(omega/c0, 2.0);
double l_out = std::pow(omega/c0, 2.0);
for (int i : cyl_dom.interior()){
op_inner.lap(i) + l_in*op_inner.value(i) = 0.0;
}
for (int i : outer_inter){
op_outer.lap(i) + l_out*op_outer.value(i) = 0.0;
}
for (int c = 0; c < interface_c_idx.size(); ++c){
int i = interface_c_idx[c]; // cyl index
int j = interface_o_idx[c]; // outer index
Vec2d pos = cyl_dom.pos(i);
double x = pos[0];
double y = pos[1];

// calculate normal derivative of the source funciton at the i-th node
// get normal, must point out of the cylinder
Vec2d normal = cyl_dom.normal(i);
std::complex<double> incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dx =
1.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dy =
1.0_i*k0*std::sin(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> dui_dn = normal[0]*din_dx + normal[1]*din_dy;

// continuity of the fields, and the incident field
op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident;

// continuity of derivatives
op_outer.neumann(j, outer_dom.normal(j))
+ (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter)
= dui_dn; // todo, change this from Neumanns to derivatives!
}
// SC-PML - where either x or y is constant
for (int i : PML){
1.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) +
((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0)
+(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0;
}
for (int i: outer_bnd){
op_outer.value(i) = 0.0;
}

</syntaxhighlight>

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics<ref name="analytical solution">M. Kerker , The Scattering of Light and Other Electromagnetic Radiation, 1969.</ref>. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{k_{0}}{k_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]
=References=

Scattering from an infinite cylinder

2026-01-22T16:24:10Z

AndrejKP: Fix the boundary condition arising from H fields (no epsilon there) and the incident wave form. Add solution reference.

Back to [[Computational electromagnetics]].

==Case==
<math>
\newcommand{\eps}{\varepsilon}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\def\doubleunderline#1{\underline{\underline{#1}}}</math>
As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
\begin{align}
\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s = 0 \qquad & \text{outside} \quad D
\label{eq:outer}
\end{align}
The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
\begin{equation}
E_z^{int} - E_z^{s} =E_z^{inc} \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n} \qquad \text{on} \quad \partial D.
\end{equation}
The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.

<syntaxhighlight lang="c++" line>
// Construct system
double l_in = epsr*std::pow(omega/c0, 2.0);
double l_out = std::pow(omega/c0, 2.0);
for (int i : cyl_dom.interior()){
op_inner.lap(i) + l_in*op_inner.value(i) = 0.0;
}
for (int i : outer_inter){
op_outer.lap(i) + l_out*op_outer.value(i) = 0.0;
}
for (int c = 0; c < interface_c_idx.size(); ++c){
int i = interface_c_idx[c]; // cyl index
int j = interface_o_idx[c]; // outer index
Vec2d pos = cyl_dom.pos(i);
double x = pos[0];
double y = pos[1];

// calculate normal derivative of the source funciton at the i-th node
// get normal, must point out of the cylinder
Vec2d normal = cyl_dom.normal(i);
std::complex<double> incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dx =
1.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> din_dy =
1.0_i*k0*std::sin(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
std::complex<double> dui_dn = normal[0]*din_dx + normal[1]*din_dy;

// continuity of the fields, and the incident field
op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident;

// continuity of derivatives
op_outer.neumann(j, outer_dom.normal(j))
+ (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter)
= dui_dn; // todo, change this from Neumanns to derivatives!
}
// SC-PML - where either x or y is constant
for (int i : PML){
1.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) +
((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0)
+(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0;
}
for (int i: outer_bnd){
op_outer.value(i) = 0.0;
}

</syntaxhighlight>

==Analytical solution==
The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
\begin{equation}
\mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
\end{equation}
The scattered field is a sum of Hankel functions so it satisfies the Sommerfeld radiation condition
\begin{equation}
\mathbf{E}^{s}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} a_{n} H_{n}^{(2)}\left(k_{0} \rho\right) e^{i n \phi}, \quad \rho \geq a
\end{equation}
and finally the total field is again a sum ob Bessel functions of the first kind
\begin{equation}
\mathbf{E}^{t}=\widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} b_{n} J_{n}\left(k_{1} \rho\right) e^{i n \phi}, \quad \rho \leq a.
\end{equation}
The boundary conditions for the field and its derivative on the boundary between the cylinder and free space surrounding it give two expressions
\begin{equation}
J_{n}\left(k_{0} a\right)+a_{n} H_{n}^{(2)}\left(k_{0} a\right)=b_{n} J_{n}\left(k_{1} a\right)
\end{equation}
and
\begin{equation}
k_0 J^{\prime}_{n}\left(k_{0} a\right)+ k_0 a_{n} H_{n}^{\prime(2)}\left(k_{0} a\right)=k_1 b_{n} J^{\prime}_{n}\left(k_{1} a\right).
\end{equation}
Which give the coefficients $a_n$ and $b_n$ in the above expansions as
\begin{equation}
a_{n}=-\frac{J_{n}\left(k_{0} a\right)-\Gamma_{h} J_{n}^{\prime}\left(k_{0} a\right)}{H_{n}^{(2)}\left(k_{0} a\right)-\Gamma_{h} H_{n}^{(2)^{\prime}}\left(k_{0} a\right)},
\end{equation}
where
\begin{equation}
\Gamma_{h}=\frac{k_{0}}{k_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
\end{equation}
and
\begin{equation}
b_n = \frac{J_n(k_0a)}{J_n(k_1a)} + a_n \frac{H_n^{(2)}(k_0a)}{J_n(k_1a)}.
\end{equation}

==Results==
[[File:FDFD_2D_dielcyl_ReEn_1.png|400px]][[File:FDFD_2D_dielcyl_ImEn_1.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_2.png|400px]][[File:FDFD_2D_dielcyl_ImEn_2.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_3.png|400px]][[File:FDFD_2D_dielcyl_ImEn_3.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_4.png|400px]][[File:FDFD_2D_dielcyl_ImEn_4.png|400px]]

[[File:FDFD_2D_dielcyl_ReEn_5.png|400px]][[File:FDFD_2D_dielcyl_ImEn_5.png|400px]]

[[File:FDFD_2D_benchmark_times_basic.png|800px]]

[[File:FDFD_2D_dx_conv2.png|800px]]

Scattering from an infinite cylinder

2026-01-22T16:22:54Z

AndrejKP:

How to build

2025-07-19T16:59:20Z

AndrejKP: Installation instructions now direct the user to master branch

=Installation=
To make this work from plain Ubuntu installation, run
<syntaxhighlight lang="bash">
sudo apt-get install git g++ python3 cmake libhdf5-serial-dev doxygen graphviz
git clone https://gitlab.com/e62Lab/medusa.git --branch master --single-branch
cd medusa
./run_tests.py -t
</syntaxhighlight>
which installs dependencies, clones the repository, goes into the root folder of
the repository and runs tests. This will build and run all tests. If this
works, you are ready to go! Otherwise install any missing packages and if it
still fails, raise an issue!
Note: If you are told the packages cannot be located, try doing a sudo apt-get update.

For instructions on how to use this library in you project, see
[[Including this library in your project]].

=Building=

List of dependencies:

* Build tools, like <code> cmake >= 2.8.12</code>, <code> g++ >= 4.8</code>, <code>make</code>, <code>python3</code>
* [https://www.hdfgroup.org/ HDF5 library] for IO
* <code>doxygen >= 1.8.8 </code> and Graphviz for generating the documentation

Out of source builds are preferred. Run
<syntaxhighlight lang="bash">
mkdir -p build
cd build
cmake ..
make
</syntaxhighlight>
Note that you only have to run <code>cmake</code> once, after that only <code>make</code> is sufficient.

Binaries are placed into <code>bin/</code> folder. Tests can be run all at once via <code>make medusa_run_tests</code>
or individually via e. g. <code>make operators_run_tests</code>.

== Linker errors ==

When trying out different classes, you might come along linker errors such as

<code>
Scanning dependencies of target cantilever_beam
[100%] Building CXX object examples/linear_elasticity/CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o
[100%] Linking CXX executable ../../../examples/linear_elasticity/cantilever_beam
/usr/bin/ld: CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o: in function `main':
cantilever_beam.cpp:(.text.startup+0x162): undefined reference to `void mm::FindBalancedSupport::operator()<mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> > >(mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&) const'
collect2: error: ld returned 1 exit status
</code>

This is expected and is the result of some optimizations of compilation time. The medusa library can actually be included in two ways: as
<code>#include <medusa/Medusa_fwd.hpp></code> or <code>#include <medusa/Medusa.hpp></code>. The first version contains the declarations of all symbols, but not all the definitions. Some of the more commonly used template instantiations are included, but by far not all. Using a template instantiation that is not precompiled will cause your program to compile fine, but will fail to link, due to the missing definitions. In this case you have a few options: include the full Medusa library (the header without the <code>_fwd</code>) and it should just work, but you will have to wait a bit longer for it to compile. Include only the missing header (in the case above <code>medusa/bits/domains/FindBalancedSupport.hpp</code>) and pay for whay you use. Or, add your instantiation among the already precompiled instantiations (located in <code>.cpp</code> files, such as e.g. [https://gitlab.com/e62Lab/medusa/-/blob/dev/src/domains/DomainDiscretization.cpp this one]).

== Building on macOS ==
This method was tested on macOS Mojave 10.14.2.

First install Xcode via App Store and then Xcode Command Line Tools with
<syntaxhighlight lang="bash">
xcode-select --install
</syntaxhighlight>

After that, install all dependencies from homebrew
<syntaxhighlight lang="bash">
brew install cmake hdf5 boost python doxygen graphviz
</syntaxhighlight>

Now you can clone and build the project using the following commands
<syntaxhighlight lang="bash">
git clone https://gitlab.com/e62Lab/medusa.git
cd medusa
mkdir build
cd build
cmake ..
cd ..
python3 run_tests.py -t
</syntaxhighlight>

This will also run all tests. If it works, you are ready to go! Otherwise install any missing packages and if it still fails, raise an issue!

==HDF5==
In order to use HDF5 IO you need the [https://www.hdfgroup.org/ HDF5 library].
You can install it easily using the command <code>sudo apt-get install libhdf5-dev</code>
or <code>sudo pacman -S hdf5</code>.

Ubuntu places (at least on older versions) hdf5 headers and libraries in a weird folder
<code>/usr/{lib, include}/x86_64-linux-gnu/hdf5/serial/</code>.

If you get an error like <code>HDF5 sample failed to compile. See errors above.</code> during <code>cmake</code> execution
then the sample hdf test file located in <code>test/test_hdf_compile.cpp</code> failed to compile. Perhaps it is good to make this file compile first,
before tackling the whole project.

If you get an error similar to <code>fatal error: hdf5.h: No such file or directory</code>,
then your compiler cannot see the HDF5 header files. Some distributions, notably (older) Ubuntu, place them into nonstandard folders
<code>/usr/include/hdf5/serial/</code> or <code>/usr/include/x86_64-linux-gnu/hdf5/serial/</code>.
Check these two folders or check your distributions hdf package for the locations of these files.
After you have determined the location, add that directory to the include directories,
using -I flag or in <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
include_directories(/usr/include/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
If you wish to fix this problem permanently, you can create soft links to the headers in your <code>/usr/include</code> directory,
by typing
<syntaxhighlight lang="bash">
sudo ln -s /usr/include/hdf5/serial/* /usr/include
</syntaxhighlight>
After this, there should be no compile time errors. If there are, please raise an issue.

If you get error similar to <code>-lhdf5 not found</code> and you have hdf5 installed,
you might have to link the libraries into a discoverable place, like <code>/usr/lib/</code>
or add the above directory to the linker path. Similarly to above, check the <code>/usr/lib/x86_64-linux-gnu/hdf5/serial/</code>
directory and look for file <code>libhdf5.a</code>. When found,
specify the location using -L flag or <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
link_directories(/usr/lib/x86_64-linux-gnu/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
or fix the problem permanently by soft-linking
<syntaxhighlight lang="bash">
sudo ln -s /usr/lib/x86_64-linux-gnu/hdf5/serial/* /usr/lib
</syntaxhighlight>

== OpenMP ==
Sometimes, OpenMP cmake errors might occure. This happens mainly due to limited multi-thread support. One can fix such issues, by adding the following code into their CMakeLists.txt

<syntaxhighlight lang="cmake">
find_package(OpenMP)
if (OPENMP_FOUND)
set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
endif()
</syntaxhighlight>

==Linear Algebra==
We use [http://eigen.tuxfamily.org/ Eigen] as our matrix library. See
[http://eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html here] for use
reference and documentation. For a quick transition from Matlab see
[http://eigen.tuxfamily.org/dox/AsciiQuickReference.txt here].

== Using Intel Math Kernel Library (MKL) ==
Install [https://software.intel.com/en-us/mkl Intel MKL] and take not of installation directories.

Proper include and link directories need to be set to use the Intel MKL.
<syntaxhighlight lang="cmake">
include_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/include) # change these to your installation path
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/lib/intel64)
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/lib/intel64)
</syntaxhighlight>

Eigen has great support for MKL. You can see more detailed instructions [https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html on their website].
To use MKL for math operations, define <code>EIGEN_USE_MKL_VML</code> when compiling. You must also link
the appropriate libraries and define <code>MKL_LP64</code> for use on a 64bit system.
With parallelism enabled, the configuration looks as follows

<syntaxhighlight lang="cmake">
target_compile_options(my_target PRIVATE "-fopenmp")
target_compile_definitions(my_target PUBLIC EIGEN_USE_MKL_VML MKL_LP64)
target_link_libraries(my_target medusa mkl_intel_lp64 mkl_intel_thread mkl_core pthread iomp5)
</syntaxhighlight>

If you have Intel Parallel Studio installed this also enables you to use the Pardiso paralle direct sparse solver through its [https://eigen.tuxfamily.org/dox/group__PardisoSupport__Module.html Eigen interface].

== Using Intel C/C++ Compiler ==

In order to use Intel's compiler when compiling Medusa, you have several standard optionas for <code>cmake</code>.
Make sure compilers and installed and in your <code>PATH</code> by running <code>which icc</code>, which
should return something like <code>/opt/intel/bin/icc</code>.

You can define the compiler when *first* calling cmake like so
<syntaxhighlight lang="bash">
cmake .. -DCMAKE_C_COMPILER=$(which icc) -DCMAKE_CXX_COMPILER=$(which icpc)
</syntaxhighlight>
If this is not your first call, remove the <code>build</code> directory and start anew.

You can also set the <code>CXX</code> and <code>CC</code> bash variables. Before calling
<code> cmake </code> for the first time you have to export the following:
<syntaxhighlight lang="bash">
export CXX="icpc"
export CC="icc"
</syntaxhighlight>

How to build

2023-06-13T11:16:47Z

AndrejKP: Temporarily changed the branch in installation instructions to "dev", until we sort out the master branch again.

=Installation=
To make this work from plain Ubuntu installation, run
<syntaxhighlight lang="bash">
sudo apt-get install git g++ python3 cmake libhdf5-serial-dev doxygen graphviz
git clone https://gitlab.com/e62Lab/medusa.git --branch dev --single-branch
cd medusa
./run_tests.py -t
</syntaxhighlight>
which installs dependencies, clones the repository, goes into the root folder of
the repository and runs tests. This will build and run all tests. If this
works, you are ready to go! Otherwise install any missing packages and if it
still fails, raise an issue!
Note: If you are told the packages cannot be located, try doing a sudo apt-get update.

For instructions on how to use this library in you project, see
[[Including this library in your project]].

=Building=

List of dependencies:

* Build tools, like <code> cmake >= 2.8.12</code>, <code> g++ >= 4.8</code>, <code>make</code>, <code>python3</code>
* [https://www.hdfgroup.org/ HDF5 library] for IO
* <code>doxygen >= 1.8.8 </code> and Graphviz for generating the documentation

Out of source builds are preferred. Run
<syntaxhighlight lang="bash">
mkdir -p build
cd build
cmake ..
make
</syntaxhighlight>
Note that you only have to run <code>cmake</code> once, after that only <code>make</code> is sufficient.

Binaries are placed into <code>bin/</code> folder. Tests can be run all at once via <code>make medusa_run_tests</code>
or individually via e. g. <code>make operators_run_tests</code>.

== Linker errors ==

When trying out different classes, you might come along linker errors such as

<code>
Scanning dependencies of target cantilever_beam
[100%] Building CXX object examples/linear_elasticity/CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o
[100%] Linking CXX executable ../../../examples/linear_elasticity/cantilever_beam
/usr/bin/ld: CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o: in function `main':
cantilever_beam.cpp:(.text.startup+0x162): undefined reference to `void mm::FindBalancedSupport::operator()<mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> > >(mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&) const'
collect2: error: ld returned 1 exit status
</code>

This is expected and is the result of some optimizations of compilation time. The medusa library can actually be included in two ways: as
<code>#include <medusa/Medusa_fwd.hpp></code> or <code>#include <medusa/Medusa.hpp></code>. The first version contains the declarations of all symbols, but not all the definitions. Some of the more commonly used template instantiations are included, but by far not all. Using a template instantiation that is not precompiled will cause your program to compile fine, but will fail to link, due to the missing definitions. In this case you have a few options: include the full Medusa library (the header without the <code>_fwd</code>) and it should just work, but you will have to wait a bit longer for it to compile. Include only the missing header (in the case above <code>medusa/bits/domains/FindBalancedSupport.hpp</code>) and pay for whay you use. Or, add your instantiation among the already precompiled instantiations (located in <code>.cpp</code> files, such as e.g. [https://gitlab.com/e62Lab/medusa/-/blob/dev/src/domains/DomainDiscretization.cpp this one]).

== Building on macOS ==
This method was tested on macOS Mojave 10.14.2.

First install Xcode via App Store and then Xcode Command Line Tools with
<syntaxhighlight lang="bash">
xcode-select --install
</syntaxhighlight>

After that, install all dependencies from homebrew
<syntaxhighlight lang="bash">
brew install cmake hdf5 boost python doxygen graphviz
</syntaxhighlight>

Now you can clone and build the project using the following commands
<syntaxhighlight lang="bash">
git clone https://gitlab.com/e62Lab/medusa.git
cd medusa
mkdir build
cd build
cmake ..
cd ..
python3 run_tests.py -t
</syntaxhighlight>

This will also run all tests. If it works, you are ready to go! Otherwise install any missing packages and if it still fails, raise an issue!

==HDF5==
In order to use HDF5 IO you need the [https://www.hdfgroup.org/ HDF5 library].
You can install it easily using the command <code>sudo apt-get install libhdf5-dev</code>
or <code>sudo pacman -S hdf5</code>.

Ubuntu places (at least on older versions) hdf5 headers and libraries in a weird folder
<code>/usr/{lib, include}/x86_64-linux-gnu/hdf5/serial/</code>.

If you get an error like <code>HDF5 sample failed to compile. See errors above.</code> during <code>cmake</code> execution
then the sample hdf test file located in <code>test/test_hdf_compile.cpp</code> failed to compile. Perhaps it is good to make this file compile first,
before tackling the whole project.

If you get an error similar to <code>fatal error: hdf5.h: No such file or directory</code>,
then your compiler cannot see the HDF5 header files. Some distributions, notably (older) Ubuntu, place them into nonstandard folders
<code>/usr/include/hdf5/serial/</code> or <code>/usr/include/x86_64-linux-gnu/hdf5/serial/</code>.
Check these two folders or check your distributions hdf package for the locations of these files.
After you have determined the location, add that directory to the include directories,
using -I flag or in <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
include_directories(/usr/include/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
If you wish to fix this problem permanently, you can create soft links to the headers in your <code>/usr/include</code> directory,
by typing
<syntaxhighlight lang="bash">
sudo ln -s /usr/include/hdf5/serial/* /usr/include
</syntaxhighlight>
After this, there should be no compile time errors. If there are, please raise an issue.

If you get error similar to <code>-lhdf5 not found</code> and you have hdf5 installed,
you might have to link the libraries into a discoverable place, like <code>/usr/lib/</code>
or add the above directory to the linker path. Similarly to above, check the <code>/usr/lib/x86_64-linux-gnu/hdf5/serial/</code>
directory and look for file <code>libhdf5.a</code>. When found,
specify the location using -L flag or <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
link_directories(/usr/lib/x86_64-linux-gnu/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
or fix the problem permanently by soft-linking
<syntaxhighlight lang="bash">
sudo ln -s /usr/lib/x86_64-linux-gnu/hdf5/serial/* /usr/lib
</syntaxhighlight>

==Linear Algebra==
We use [http://eigen.tuxfamily.org/ Eigen] as our matrix library. See
[http://eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html here] for use
reference and documentation. For a quick transition from Matlab see
[http://eigen.tuxfamily.org/dox/AsciiQuickReference.txt here].

== Using Intel Math Kernel Library (MKL) ==
Install [https://software.intel.com/en-us/mkl Intel MKL] and take not of installation directories.

Proper include and link directories need to be set to use the Intel MKL.
<syntaxhighlight lang="cmake">
include_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/include) # change these to your installation path
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/lib/intel64)
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/lib/intel64)
</syntaxhighlight>

Eigen has great support for MKL. You can see more detailed instructions [https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html on their website].
To use MKL for math operations, define <code>EIGEN_USE_MKL_VML</code> when compiling. You must also link
the appropriate libraries and define <code>MKL_LP64</code> for use on a 64bit system.
With parallelism enabled, the configuration looks as follows

<syntaxhighlight lang="cmake">
target_compile_options(my_target PRIVATE "-fopenmp")
target_compile_definitions(my_target PUBLIC EIGEN_USE_MKL_VML MKL_LP64)
target_link_libraries(my_target medusa mkl_intel_lp64 mkl_intel_thread mkl_core pthread iomp5)
</syntaxhighlight>

If you have Intel Parallel Studio installed this also enables you to use the Pardiso paralle direct sparse solver through its [https://eigen.tuxfamily.org/dox/group__PardisoSupport__Module.html Eigen interface].

== Using Intel C/C++ Compiler ==

In order to use Intel's compiler when compiling Medusa, you have several standard optionas for <code>cmake</code>.
Make sure compilers and installed and in your <code>PATH</code> by running <code>which icc</code>, which
should return something like <code>/opt/intel/bin/icc</code>.

You can define the compiler when *first* calling cmake like so
<syntaxhighlight lang="bash">
cmake .. -DCMAKE_C_COMPILER=$(which icc) -DCMAKE_CXX_COMPILER=$(which icpc)
</syntaxhighlight>
If this is not your first call, remove the <code>build</code> directory and start anew.

You can also set the <code>CXX</code> and <code>CC</code> bash variables. Before calling
<code> cmake </code> for the first time you have to export the following:
<syntaxhighlight lang="bash">
export CXX="icpc"
export CC="icc"
</syntaxhighlight>

Computation of shape functions

2022-06-16T12:58:39Z

AndrejKP: Typos

Shape functions, often called stencil weights, are functionals, approximating linear partial differential operators at a chosen point.
A shape function $\b{\chi}_{\mathcal{L}} (\vec{p})$ for a linear partial differential operator $\mathcal{L}$
in a point $p$ is a vector of stencil weights, such that
\[
(\mathcal{L} u)(\vec{p}) \approx \b{\chi}_{\mathcal{L}}(\vec{p})^\T \b{u},
\]
where $\b{u}$ is a vector of function values $u(x_i)$ evaluated in the support of $\vec{p}$.
For usage examples see the section on [[Meshless Local Strong Form Method (MLSM)]] or some [[ Poisson's equation | examples ]].

For computation of shape functions augmented with monomials, see the [[ Radial basis function-generated finite differences (RBF-FD) ]].

== Derivation ==
Let $\mathcal{L}$ be a linear partial differential operator.
We wish to approximate $(\mathcal{L}u)(\vec{p})$ using only values of $u$ in support of $p$. To do so, we construct a
WLS approximation (see [[Weighted Least Squares (WLS)]]) and compute the coefficients $\b{\alpha}$, such that $\hat{u}(\vec{p}) = \b{b}(\vec{p})^\T \b{\alpha}$.
The coefficients are computed by solving a least-squares system $WB \b{\alpha} = W\b{u}$, (see [[Solving overdetermined systems]])
writing
\[ \b{\alpha} = (B^\T W^2 B^\T)^{-1} B^\T W^2 \b{u}. \]
Writing down the approximation function $\hat{u}$ we get
\[
\hat u (\vec{p}) = \b{b}(\vec{p})^\T \b{\alpha} = \b{b}(\vec{p})^\T (WB)^+W\b{u}
\]
and applying operator $\mathcal{L}$ we get
\[
\mathcal{L} \hat u (\vec{p}) = (\mathcal{L}\b{b})(\vec{p})^\T \b{\alpha} =
(\mathcal{L}\b{b})(\vec{p})^\T (B^\T W^2 B)^{-1} B^\T W^2 \b{u}
= \b{\chi}_{\mathcal{L}}(\vec{p})^\T \b{u}
\]

Depending on what data we have available, we can either compute $\b \alpha$ and obtain a continuous approximation $\hat{u}$,
to which $\mathcal{L}$ can be applied at any point, or, if function values $\b u$ are unknown, we can choose to compute
$\b \chi_{\mathcal{L}}$ at a point $\vec p$ and approximate the operator for ''any'' set of function values.
the operator $\mathcal{L}$ at $\vec p$ is approximated as
\[
(\mathcal{L}u)(\vec{p}) \approx (\mathcal{L} \hat{u})(p) = \b{\chi}_{\mathcal{L}}(\vec{p})^\T \b{u}.
\]
Vector $\b{\chi}$ is a vector of stencil weights, also called a ''shape function''. The name comes historically from FEM, and it incorporates the idea of being able to take all the information
about shape of the local domain and choice of approximation and store it in a single vector, being able to approximate
a function value from given support values $\b{u}$ with a single dot product. For any values $\b{u}$, value $\b{\chi}(p) \b{u}$
gives us the approximation of $(\mathcal{L}u)(\vec{p})$.
Mathematically speaking, $\b{\chi}(\vec{p})$ is a Riesz representation of a functional, $\b{\chi}(\vec{p})\colon \R^n \to \R$, mapping $n$-tuples of known function values to
approximations of $\mathcal{L}$ in point $\vec{p}$.

For example, take a $1$-dimensional case for approximation of derivatives with weight equal to $1$ and $n=m=3$, with equally spaced support values at distances $h$.
We wish to approximate $u''$ in the middle support point, just by making a weighted sum of the values, something like the finite difference
\[ u'' \approx \frac{u_1 - 2u_2 + u_3}{h^2}. \]
This is exactly the same formula as we would have come to by computing $\b{\chi}$, except that our approach is a lot more general. But one should think about
$\b{\chi}$ as one would about the finite difference scheme, it is a rule, telling us how to compute the derivative.
\[ u''(s_2) \approx \underbrace{\begin{bmatrix} \frac{1}{h^2} & \frac{-2}{h^2} & \frac{1}{h^2} \end{bmatrix}}_{\b{\chi}^\T} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} \]

The fact that $\b{\chi}$ is independent of the function values $\b{u}$ but depends only on domain geometry means that
'''we can just compute the shape functions $\b{\chi}$ for points of interest and then approximate any linear operator
of any function, given its values, very fast, using only a single dot product.'''

== Special case when $B$ is square and invertible ==
The expression
\[ \b{\chi}_{\mathcal{L}}(\vec{p})^\T = (\mathcal{L}\b{b})(\vec{p})^\T (B^\T W^2 B)^{-1} B^\T W^2 \]
can be simplified to
\[
\b \chi_{\mathcal{L}}(\vec{p})^\T = (\mathcal{L}\b{b})(\vec{p})^\T B^{-1}.
\]

This gives the approximation as $(\mathcal{L} \hat{u})(p) = (\mathcal{L}\b{b})(\vec{p})^\T B^{-1} \b u$.

Coefficients $\b \alpha$ are the solutions of $B\b \alpha = \b u$. Shape function $\b{\chi}_{\mathcal{L}}(\vec{p})$ is defined as
$$
\b{\chi}_{\mathcal{L}}(\vec{p})^\T = (\mathcal{L}\b{b})(\vec{p})^\T B^{-1},
$$
which can be rewritten as a solution of
$$ B^\T \b{\chi}_{\mathcal{L}}(\vec{p}) = (\mathcal{L}\b{b})(\vec{p}). $$

This formulation is more suitable for numerical solving and also has a deeper interpretation, as presented in the next section.

== Another view on derivation of shape functions ==
Solving the system $B^\T \b{\chi}_{\mathcal{L}} (\vec{p})^\T = (\mathcal{L}\b{b})(\vec p)$ can be interpreted as using the method of undetermined coefficients.

We wish to find such a combination of function values evaluated in support points $\vec s_i$ that approximates the value of an operator $(\mathcal{L}u)(\vec p)$, i.e., we are looking for coefficients $\chi_{\mathcal{L}, i}$, such that
\[ \sum_{i=1}^n \chi_{\mathcal{L}, i} u_i \approx (\mathcal{L}u)(p). \]

We require that $\chi_{\mathcal{L}, i}$ be such, that above approximations holds for a certain space of functions, spanned by the basis $\{b_j\}_{j=1}^m$.
By choosing the basis functions $b_j$ we decide for which functions the aproximation $\b{\chi}_{\mathcal{L}, \vec{p}} \b{u} \approx (\mathcal{L}u)(\vec p)$ should be exact and we write the equations
$\sum_{i=1}^n \chi_{\mathcal{L}, i} b_j(\vec s_i) = (\mathcal{L}b_j)(\vec p)$ in a system:
\[
\begin{bmatrix}
b_1(\vec s_1) & \cdots & b_1(\vec s_n) \\
\vdots & \ddots & \vdots \\
b_m(\vec s_1) & \cdots & b_m(\vec s_n) \\
\end{bmatrix}
\begin{bmatrix}
\chi_{\mathcal{L}, 1} \\ \vdots \\ \chi_{\mathcal{L}, n}
\end{bmatrix}
=
\begin{bmatrix}
(\mathcal{L} b_1)(\vec{p}) \\ \vdots \\ (\mathcal{L} b_m)(\vec{p})
\end{bmatrix}
\]
Note that the matrix above is precisely the transpose of the matrix $B$ from the [[Weighted Least Squares (WLS) ]] approximation. In this sense, the computation of
shape functions is dual to computation of basis coefficients. The system is written in matrix form as $B^\T \b{\chi}_{\mathcal{L}} (\vec{p})^\T = (\mathcal{L}\b{b})(\vec p)$.

For $n=m$ case as in the previous section, the solution is the sought shape function $\b{\chi}_{\mathcal L}$.

If we do not specify enough functions $b_j$, the system above is underdetermined.
A common solution is to minimize the weighted norm of $\b \chi_{\mathcal L}$.
This gives us a minimization problem
$$
\min_{\b \chi_{\mathcal L}} \sum_{i=1}^n (\chi_{\mathcal{L}, i}/w_i)^2, \text{ such that } B^\T \b{\chi}_{\mathcal{L}} (\vec{p})^\T = (\mathcal{L}\b{b})(\vec p).
$$
If the central weights are bigger, this also forces the central coefficients to be larger.

Using the solution derived in [[Solving underdetermined systems]] and substituting inverse weights $W = W^{-1}$ and matrix $A = B^\T$, we get the solution
$$
\b{\chi}_{\mathcal{L}} (\vec{p}) = (W^{-1})^{-2} (B^\T)^\T (B^\T (W^{-1})^{-2}(B^\T)^\T)^{-1} (\mathcal{L}\b{b})(\vec p) = W^2 B (B^\T W^2 B)^{-1} (\mathcal{L}\b{b})(\vec p).
$$
After transposing, we get
$$
\b{\chi}_{\mathcal{L}} (\vec{p})^\T = (\mathcal{L}\b{b})(\vec p)^\T (B^\T W^2 B)^{-1} B^\T W^2,
$$
which is exactly the same formula as derived using [[ Weighted Least Squares (WLS) ]] approximation.

Just as we get the same stencil weights using interpolation and the method of undetermined coefficients, we get the same result using
WLS approximation or weighted norm minimization of shape functions.

== Numerical calculation of the shape functions ==
Numerically cheaper and more stable way is to translate the problem of inverting the matrix to solving a linear system of equations.
The alternative view above is useful in this aspect.

=== Invertible $B$ ===

If $B$ is invertible, then $\b{\chi}_{\mathcal{L}}(\vec{p})$ can be thought as a solution of a system $B^\T\b{\chi}_{\mathcal{L}}(\vec{p}) = (\mathcal{L}\b{b})(\vec p)$, which can be solved using LU or Cholesky decomposition for example.
For more on this, see [[Solving linear systems]].

=== General case ===
In general, the system $B^\T\b{\chi}_{\mathcal{L}}(\vec{p}) = (\mathcal{L}\b{b})(\vec p)$ is underdetermined and must be solved in the least-weighted-norm sense.
This means solving the underdetermined system $B^\T W y = (\mathcal{L}\b{b})(\vec p)$ and computing $\b{\chi}_{\mathcal{L}}(\vec{p}) = Wy$.
For more details on this, see [[Solving underdetermined systems]].

Depending on our knowledge of the
problem, we can use Cholesky ($B^\T W$ is positive definite), $LDL^\T$ if it is symmetric, $LU$ for a general square matrix, $QR$ for full rank overdetermined system and SVD for a general system.

If more shapes need to be calculated using the same matrix $B^\T W$ and only different right hand sides (as is often the case, if shapes for more than one operator are computed),
it can be done efficiently by storing the decomposition of $B^\T W$.

== Linearity and reduction to only computing the shapes of the operator space basis ==
Observe that the expression for the shape functions
\[ \b{\chi}_{\mathcal{L}}(\vec{p}) = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+W \]
is linear in $\mathcal{L}$. This follows from the fact that if $\mathcal{L} = \mathcal{L}_1 + \mathcal{L}_2$, the vector of
derivative values $(\mathcal{L}\b{b})(\vec{p}) = ((\mathcal{L}_1 + \mathcal{L}_2)\b{b})(\vec{p}) =
(\mathcal{L}_1\b{b} + \mathcal{L}_2\b{b})(\vec{p}) = (\mathcal{L}_1\b{b})(\vec p) + (\mathcal{L}_2\b{b})(\vec p)$
and consequently
$$\b{\chi}_{\mathcal{L}}(\vec{p}) = \b{\chi}_{\mathcal{L}_1}(\vec{p}) + \b{\chi}_{\mathcal{L}_2}(\vec{p}).$$

Say that your equation constitutes of operators up to a certain $k$-th order.
Then it is sufficient to compute only shape functions for the basis of that operator space, namely
\[
\{\b{\chi}_{\frac{\partial}{\partial x^\alpha}, p}, 1 \leq |\alpha| \leq k \}, \text{ where } \frac{\partial}{\partial x^\alpha} = \frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1}\cdots \partial x_d^{\alpha_d}} \text{ and } |\alpha| = \sum_{i=1}^d \alpha_i.
\]
Having the shape functions $\b{\chi}_{\frac{\partial}{\partial x^\alpha}, p}$ available, we write $\mathcal{L}$ in the above basis in point $p$, $\mathcal{L} = \displaystyle \sum_{1 \leq |\alpha| \leq k} a_\alpha(p) \frac{\partial}{\partial x^\alpha}$
and evaluate it's shape function as
\[
\b{\chi}_{\mathcal{L}, \vec{p}} = \sum_{1 \leq |\alpha| \leq k} a_\alpha(p) \b{\chi}_{\frac{\partial}{\partial x^\alpha}, p}.
\]
This saves us space and gives certain elegance to the implementation.

Weighted Least Squares (WLS)

2022-06-16T12:57:59Z

AndrejKP: Typos

One of the most important building blocks of the meshless methods is the Weighted/Moving Least Squares (WLS/MLS) approximation , which is implemented in the <code>WLS</code> class.
It is also often called GFDM (generalized finite differences method).

= Notation cheat sheet =
\begin{align*}
m \in \N & \dots \text{number of basis functions} \\
n \geq m \in \N & \dots \text{number of points in support domain} \\
k \in \mathbb{N} & \dots \text{dimensionality of vector space} \\
\vec s_j \in \R^k & \dots \text{point in support domain } \quad j=1,\dots,n \\
u_j \in \R & \dots \text{value of function to approximate in }\vec{s}_j \quad j=1,\dots,n \\
\vec p \in \R^k & \dots \text{center point of approximation} \\
b_i\colon \R^k \to \R & \dots \text{basis functions } \quad i=1,\dots,m \\
B_{j, i} \in \R & \dots \text{value of basis functions in support points } b_i(\vec{s}_j) \quad j=1,\dots,n, \quad i=1,\dots,m\\
\omega \colon \R^k \to \R & \dots \text{weight function} \\
w_j \in \R & \dots \text{weights } \omega(\vec{s}_j-\vec{p}) \quad j=1,\dots,n \\
\alpha_i \in \R & \dots \text{expansion coefficients around point } \vec{p} \quad i=1,\dots,m \\
\hat u\colon \R^k \to \R & \dots \text{approximation function (best fit)} \\
\chi_j \in \R & \dots \text{shape coefficient for point }\vec{p} \quad j=1,\dots,n \\
\end{align*}

We will also use $\b{s}, \b{u}, \b{b}, \b{\alpha}, \b{\chi} $ to denote a column of corresponding values,
$W$ as a $n\times n$ diagonal matrix filled with $w_j$ on the diagonal and $B$ as a $n\times m$ matrix filled with $B_{j, i} = b_i(\vec s_j)$.

= Definition of local approximation =
<figure id="fig:1DWLS">
[[File:image_1avhdsfej1b9cao01029m1e13o69.png|600px|thumb|upright=2|alt=1D MLS example|<caption>Example of 1D WLS approximation </caption>]]
</figure>
Our wish is to approximate an unknown function $u\colon \R^k \to \R$ while knowing $n$ values $u(\vec{s}_j) := u_j$.
The vector of known values will be denoted by $\b{u}$ and the vector of coordinates where those values were achieved by $\b{s}$.
Note that $\b{s}$ is not a vector in the usual sense since its components $\vec{s}_j$ are elements of $\R^k$, but we will call it vector anyway.
The values of $\b{s}$ are called ''nodes'' or ''support nodes'' or ''support''. The known values $\b{u}$ are also called ''support values''.

In general, an approximation function around point $\vec{p}\in\R^k$ can be
written as \[\hat{u} (\vec{x}) = \sum_{i=1}^m \alpha_i b_i(\vec{x}) = \b{b}(\vec{x})^\T \b{\alpha} \]
where $\b{b} = (b_i)_{i=1}^m$ is a set of ''basis functions'', $b_i\colon \R^k \to\R$, and $\b{\alpha} = (\alpha_i)_{i=1}^m$ are the unknown coefficients.

In MLS the goal is to minimize the error of approximation in given values, $\b{e} = \hat u(\b{s}) - \b{u}$
between the approximation function and target function in the known points $\b{x}$. The error can also be written as $B\b{\alpha} - \b{u}$,
where $B$ is rectangular matrix of dimensions $n \times m$ with rows containing basis function evaluated in points $\vec{s}_j$.
\[ B =
\begin{bmatrix}
b_1(\vec{s}_1) & \ldots & b_m(\vec{s}_1) \\
\vdots & \ddots & \vdots \\
b_1(\vec{s}_n) & \ldots & b_m(\vec{s}_n)
\end{bmatrix} =
[b_i(\vec{s}_j)]_{j=1,i=1}^{n,m} = [\b{b}(\vec{s}_j)^\T]_{j=1}^n. \]

We can choose to minimize any norm of the error vector $e$
and usually choose to minimize the $2$-norm or square norm \[ \|\b{e}\| = \|\b{e}\|_2 = \sqrt{\sum_{j=1}^n e_j^2}. \]
Commonly, we also choose to minimize a weighted norm
<ref>Note that our definition is a bit unusual, usually weights are not
squared with the values. However, we do this to avoid computing square
roots when doing MLS. If you are used to the usual definition,
consider the weight to be $\omega^2$.</ref>
instead \[ \|\b{e}\|_{2,w} = \|\b{e}\|_w = \sqrt{\sum_{j=1}^n (w_j e_j)^2}. \]
The ''weights'' $w_i$ are assumed to be non negative and are assembled in a vector $\b{w}$ or a matrix $W = \operatorname{diag}(\b{w})$ and usually obtained from a weight function.
A ''weight function'' is a function $\omega\colon \R^k \to[0,\infty)$. We calculate $w_j$ as $w_i := \omega(\vec{p}-\vec{s}_j)$, so
good choices for $\omega$ are function which have higher values close to $0$ (making closer nodes more important), like the normal distribution.
If we choose $\omega \equiv 1$, we get the unweighted version.

A choice of minimizing the square norm gave this method its name - Least Squares approximation. If we use the weighted version, we get the Weighted Least Squares or WLS.
In the most general case we wish to minimize
\[ \|\b{e}\|_{2,w}^2 = \b{e}^\T W^2 \b{e} = (B\b{\alpha} - \b{u})^\T W^2(B\b{\alpha} - \b{u}) = \sum_j^n w_j^2 (\hat{u}(\vec{s}_j) - u_j)^2 \]

The problem of finding the coefficients $\b{\alpha}$ that minimize the error $\b{e}$ numerically is described in [[Solving overdetermined systems]].

Here, we simply write the solution $\b \alpha = (WB)^+ Wu = (B^\T W^2 B)^{-1} B^\T W^2 u$.

In our WLS engine we use SVD with regularization by default.

= Weighted Least Squares =
Weighted least squares approximation is the simplest version of the procedure described above. Given support $\b{s}$, values $\b{u}$
and an anchor point $\vec{p}$, we calculate the coefficients $\b{\alpha}$ using one of the above methods.
Then, to approximate a function in the neighbourhood of $\vec p$ we use the formula
\[
\hat{u}(\vec x) = \b{b}(\vec x)^\T \b{\alpha} = \sum_{i=1}^m \alpha_i b_i(\vec x).
\]

To approximate the derivative $\frac{\partial u}{\partial x_i}$, or any linear partial differential operator $\mathcal L$ on $u$, we
simply take the same linear combination of transformed basis functions $\mathcal L b_i$. We have considered coefficients $\alpha_i$ to be
constant and applied the linearity.
\[
\widehat{\mathcal L u}(\vec x) = \sum_{i=1}^m \alpha_i (\mathcal L b_i)(\vec x) = (\mathcal L \b b(\vec x))^\T \b \alpha.
\]

= Moving Least Squares =

<figure id="fig:comparisonMLSandWLS">
[[File:mlswls.svg|thumb|upright=2|<caption>Comparison of WLS and MLS approximation</caption>]]
</figure>

When using WLS the approximation gets worse as we move away from the central point $\vec{p}$.
This is partially due to not being in the center of the support any more and partially due to weight
being distributed in such a way to assign more importance to nodes closer to $\vec{p}$.

We can battle this problem in two ways: when we wish to approximate in a new point that is sufficiently far
away from $\vec{p}$ we can compute new support, recompute the new coefficients $\b{\alpha}$ and approximate again.
This is very costly and we would like to avoid that. A partial fix is to keep support the same, but only
recompute the weight vector $\b{w}$, that will now assign weight values to nodes close to the new point.
We still need to recompute the coefficients $\b{\alpha}$, however we avoid the cost of setting up new support
and function values and recomputing $B$. This approach is called Moving Least Squares due to recomputing
the weighted least squares problem whenever we move the point of approximation.

Note that if our weight is constant or if $n = m$, when approximation reduces to interpolation, the weights do not play
any role and this method is redundant. In fact, its benefits arise when supports are rather large.

See <xr id="fig:comparisonMLSandWLS"/> for comparison between MLS and WLS approximations. MLS approximation remains close to
the actual function while still inside the support domain, while WLS approximation becomes bad when
we come out of the reach of the weight function.

{{reflist}}

Weighted Least Squares (WLS)

2022-06-16T12:55:27Z

AndrejKP: b_i (s_j - p) -> b_i (s_j) in notation cheat sheet (seems like what is used from then on)

One of the most important building blocks of the meshless methods is the Weighted/Moving Least Squares (WLS/MLS) approximation , which is implemented in the <code>WLS</code> class.
It is also often called GFDM (generalized finite differences method).

= Notation cheat sheet =
\begin{align*}
m \in \N & \dots \text{number of basis functions} \\
n \geq m \in \N & \dots \text{number of points in support domain} \\
k \in \mathbb{N} & \dots \text{dimensionality of vector space} \\
\vec s_j \in \R^k & \dots \text{point in support domain } \quad j=1,\dots,n \\
u_j \in \R & \dots \text{value of function to approximate in }\vec{s}_j \quad j=1,\dots,n \\
\vec p \in \R^k & \dots \text{center point of approximation} \\
b_i\colon \R^k \to \R & \dots \text{basis functions } \quad i=1,\dots,m \\
B_{j, i} \in \R & \dots \text{value of basis functions in support points } b_i(\vec{s}_j) \quad j=1,\dots,n, \quad i=1,\dots,m\\
\omega \colon \R^k \to \R & \dots \text{weight function} \\
w_j \in \R & \dots \text{weights } \omega(\vec{s}_j-\vec{p}) \quad j=1,\dots,n \\
\alpha_i \in \R & \dots \text{expansion coefficients around point } \vec{p} \quad i=1,\dots,m \\
\hat u\colon \R^k \to \R & \dots \text{approximation function (best fit)} \\
\chi_j \in \R & \dots \text{shape coefficient for point }\vec{p} \quad j=1,\dots,n \\
\end{align*}

We will also use $\b{s}, \b{u}, \b{b}, \b{\alpha}, \b{\chi} $ to denote a column of corresponding values,
$W$ as a $n\times n$ diagonal matrix filled with $w_j$ on the diagonal and $B$ as a $n\times m$ matrix filled with $B_{j, i} = b_i(\vec s_j)$.

= Definition of local approximation =
<figure id="fig:1DWLS">
[[File:image_1avhdsfej1b9cao01029m1e13o69.png|600px|thumb|upright=2|alt=1D MLS example|<caption>Example of 1D WLS approximation </caption>]]
</figure>
Our wish is to approximate an unknown function $u\colon \R^k \to \R$ while knowing $n$ values $u(\vec{s}_j) := u_j$.
The vector of known values will be denoted by $\b{u}$ and the vector of coordinates where those values were achieved by $\b{s}$.
Note that $\b{s}$ is not a vector in the usual sense since its components $\vec{s}_j$ are elements of $\R^k$, but we will call it vector anyway.
The values of $\b{s}$ are called ''nodes'' or ''support nodes'' or ''support''. The known values $\b{u}$ are also called ''support values''.

In general, an approximation function around point $\vec{p}\in\R^k$ can be
written as \[\hat{u} (\vec{x}) = \sum_{i=1}^m \alpha_i b_i(\vec{x}) = \b{b}(\vec{x})^\T \b{\alpha} \]
where $\b{b} = (b_i)_{i=1}^m$ is a set of ''basis functions'', $b_i\colon \R^k \to\R$, and $\b{\alpha} = (\alpha_i)_{i=1}^m$ are the unknown coefficients.

In MLS the goal is to minimize the error of approximation in given values, $\b{e} = \hat u(\b{s}) - \b{u}$
between the approximation function and target function in the known points $\b{x}$. The error can also be written as $B\b{\alpha} - \b{u}$,
where $B$ is rectangular matrix of dimensions $n \times m$ with rows containing basis function evaluated in points $\vec{s}_j$.
\[ B =
\begin{bmatrix}
b_1(\vec{s}_1) & \ldots & b_m(\vec{s}_1) \\
\vdots & \ddots & \vdots \\
b_1(\vec{s}_n) & \ldots & b_m(\vec{s}_n)
\end{bmatrix} =
[b_i(\vec{s}_j)]_{j=1,i=1}^{n,m} = [\b{b}(\vec{s}_j)^\T]_{j=1}^n. \]

We can choose to minimize any norm of the error vector $e$
and usually choose to minimize the $2$-norm or square norm \[ \|\b{e}\| = \|\b{e}\|_2 = \sqrt{\sum_{j=1}^n e_j^2}. \]
Commonly, we also choose to minimize a weighted norm
<ref>Note that our definition is a bit unusual, usually weights are not
squared with the values. However, we do this to avoid computing square
roots when doing MLS. If you are used to the usual definition,
consider the weight to be $\omega^2$.</ref>
instead \[ \|\b{e}\|_{2,w} = \|\b{e}\|_w = \sqrt{\sum_{j=1}^n (w_j e_j)^2}. \]
The ''weights'' $w_i$ are assumed to be non negative and are assembled in a vector $\b{w}$ or a matrix $W = \operatorname{diag}(\b{w})$ and usually obtained from a weight function.
A ''weight function'' is a function $\omega\colon \R^k \to[0,\infty)$. We calculate $w_j$ as $w_i := \omega(\vec{p}-\vec{s}_j)$, so
good choices for $\omega$ are function which have higher values close to $0$ (making closer nodes more important), like the normal distribution.
If we choose $\omega \equiv 1$, we get the unweighted version.

A choice of minimizing the square norm gave this method its name - Least Squares approximation. If we use the weighted version, we get the Weighted Least Squares or WLS.
In the most general case we wish to minimize
\[ \|\b{e}\|_{2,w}^2 = \b{e}^\T W^2 \b{e} = (B\b{\alpha} - \b{u})^\T W^2(B\b{\alpha} - \b{u}) = \sum_j^n w_j^2 (\hat{u}(\vec{s}_j) - u_j)^2 \]

The problem of finding the coefficients $\b{\alpha}$ that minimize the error $\b{e}$ numerically is described in [[Solving overdetermined systems]].

Here, we simply write the solution $\b \alpha = (WB)^+ Wu = (B^\T W^2 B)^{-1} B^\T W^2 u$.

In our WLS engine we use SVD with regularization by default.

= Weighted Least Squares =
Weighted least squares approximation is the simplest version of the procedure described above. Given support $\b{s}$, values $\b{u}$
and an anchor point $\vec{p}$, we calculate the coefficients $\b{\alpha}$ using one of the above methods.
Then, to approximate a function in the neighbourhood of $\vec p$ we use the formula
\[
\hat{u}(\vec x) = \b{b}(\vec x)^\T \b{\alpha} = \sum_{i=1}^m \alpha_i b_i(\vec x).
\]

To approximate the derivative $\frac{\partial u}{\partial x_i}$, or any linear partial differential operator $\mathcal L$ on $u$, we
simply take the same linear combination of transformed basis functions $\mathcal L b_i$. We have considered coefficients $\alpha_i$ to be
constant and applied the linearity.
\[
\widehat{\mathcal L u}(\vec x) = \sum_{i=1}^m \alpha_i (\mathcal L b_i)(\vec x) = (\mathcal L \b b(\vec x))^\T \b \alpha.
\]

= Moving Least Squares =

<figure id="fig:comparisonMLSandWLS">
[[File:mlswls.svg|thumb|upright=2|<caption>Comparison of WLS and MLS approximation</caption>]]
</figure>

When using WLS the approximation gets worse as we move away from the central point $\vec{p}$.
This is partially due to not being in the center of the support any more and partially due to weight
being distributed in such a way to assign more importance to nodes closer to $\vec{p}$.

We can battle this problem in two ways: when we wish to approximate in a new point that is sufficiently far
away from $\vec{p}$ we can compute new support, recompute the new coefficients $\b{\alpha}$ and approximate again.
This is very costly and we would like to avoid that. A partial fix is to keep support the same, but only
recompute the weight vector $\b{w}$, that will now assign weight values to nodes close to the new point.
We still need to recompute the coefficients $\b{\alpha}$, however we avoid the cost of setting up new support
and function values and recomputing $B$. This approach is called Moving Least Squares due to recomputing
the weighted least squares problem whenever we move the point of approximation.

Note that if out weight is constant or if $n = m$, when approximation reduces to interpolation, the weights do not play
any role and this method is redundant. In fact, its benefits arise when supports are rather large.

See <xr id="fig:comparisonMLSandWLS"/> for comparison between MLS and WLS approximations. MLS approximation remains close to
actual function while still inside the support domain, while WLS approximation becomes bad when
we come out of the reach of the weight function.

{{reflist}}

How to build

2022-06-16T12:52:28Z

AndrejKP: Added the mention of the fix for "packages cannot be located" (at least it worked for me)

=Installation=
To make this work from plain Ubuntu installation, run
<syntaxhighlight lang="bash">
sudo apt-get install git g++ python3 cmake libhdf5-serial-dev doxygen graphviz
git clone https://gitlab.com/e62Lab/medusa.git --branch master --single-branch
cd medusa
./run_tests.py -t
</syntaxhighlight>
which installs dependencies, clones the repository, goes into the root folder of
the repository and runs tests. This will build and run all tests. If this
works, you are ready to go! Otherwise install any missing packages and if it
still fails, raise an issue!
Note: If you are told the packages cannot be located, try doing a sudo apt-get update.

For instructions on how to use this library in you project, see
[[Including this library in your project]].

=Building=

List of dependencies:

* Build tools, like <code> cmake >= 2.8.12</code>, <code> g++ >= 4.8</code>, <code>make</code>, <code>python3</code>
* [https://www.hdfgroup.org/ HDF5 library] for IO
* <code>doxygen >= 1.8.8 </code> and Graphviz for generating the documentation

Out of source builds are preferred. Run
<syntaxhighlight lang="bash">
mkdir -p build
cd build
cmake ..
make
</syntaxhighlight>
Note that you only have to run <code>cmake</code> once, after that only <code>make</code> is sufficient.

Binaries are placed into <code>bin/</code> folder. Tests can be run all at once via <code>make medusa_run_tests</code>
or individually via e. g. <code>make operators_run_tests</code>.

== Linker errors ==

When trying out different classes, you might come along linker errors such as

<code>
Scanning dependencies of target cantilever_beam
[100%] Building CXX object examples/linear_elasticity/CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o
[100%] Linking CXX executable ../../../examples/linear_elasticity/cantilever_beam
/usr/bin/ld: CMakeFiles/cantilever_beam.dir/cantilever_beam.cpp.o: in function `main':
cantilever_beam.cpp:(.text.startup+0x162): undefined reference to `void mm::FindBalancedSupport::operator()<mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> > >(mm::DomainDiscretization<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&) const'
collect2: error: ld returned 1 exit status
</code>

This is expected and is the result of some optimizations of compilation time. The medusa library can actually be included in two ways: as
<code>#include <medusa/Medusa_fwd.hpp></code> or <code>#include <medusa/Medusa.hpp></code>. The first version contains the declarations of all symbols, but not all the definitions. Some of the more commonly used template instantiations are included, but by far not all. Using a template instantiation that is not precompiled will cause your program to compile fine, but will fail to link, due to the missing definitions. In this case you have a few options: include the full Medusa library (the header without the <code>_fwd</code>) and it should just work, but you will have to wait a bit longer for it to compile. Include only the missing header (in the case above <code>medusa/bits/domains/FindBalancedSupport.hpp</code>) and pay for whay you use. Or, add your instantiation among the already precompiled instantiations (located in <code>.cpp</code> files, such as e.g. [https://gitlab.com/e62Lab/medusa/-/blob/dev/src/domains/DomainDiscretization.cpp this one]).

== Building on macOS ==
This method was tested on macOS Mojave 10.14.2.

First install Xcode via App Store and then Xcode Command Line Tools with
<syntaxhighlight lang="bash">
xcode-select --install
</syntaxhighlight>

After that, install all dependencies from homebrew
<syntaxhighlight lang="bash">
brew install cmake hdf5 boost python doxygen graphviz
</syntaxhighlight>

Now you can clone and build the project using the following commands
<syntaxhighlight lang="bash">
git clone https://gitlab.com/e62Lab/medusa.git
cd medusa
mkdir build
cd build
cmake ..
cd ..
python3 run_tests.py -t
</syntaxhighlight>

This will also run all tests. If it works, you are ready to go! Otherwise install any missing packages and if it still fails, raise an issue!

==HDF5==
In order to use HDF5 IO you need the [https://www.hdfgroup.org/ HDF5 library].
You can install it easily using the command <code>sudo apt-get install libhdf5-dev</code>
or <code>sudo pacman -S hdf5</code>.

Ubuntu places (at least on older versions) hdf5 headers and libraries in a weird folder
<code>/usr/{lib, include}/x86_64-linux-gnu/hdf5/serial/</code>.

If you get an error like <code>HDF5 sample failed to compile. See errors above.</code> during <code>cmake</code> execution
then the sample hdf test file located in <code>test/test_hdf_compile.cpp</code> failed to compile. Perhaps it is good to make this file compile first,
before tackling the whole project.

If you get an error similar to <code>fatal error: hdf5.h: No such file or directory</code>,
then your compiler cannot see the HDF5 header files. Some distributions, notably (older) Ubuntu, place them into nonstandard folders
<code>/usr/include/hdf5/serial/</code> or <code>/usr/include/x86_64-linux-gnu/hdf5/serial/</code>.
Check these two folders or check your distributions hdf package for the locations of these files.
After you have determined the location, add that directory to the include directories,
using -I flag or in <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
include_directories(/usr/include/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
If you wish to fix this problem permanently, you can create soft links to the headers in your <code>/usr/include</code> directory,
by typing
<syntaxhighlight lang="bash">
sudo ln -s /usr/include/hdf5/serial/* /usr/include
</syntaxhighlight>
After this, there should be no compile time errors. If there are, please raise an issue.

If you get error similar to <code>-lhdf5 not found</code> and you have hdf5 installed,
you might have to link the libraries into a discoverable place, like <code>/usr/lib/</code>
or add the above directory to the linker path. Similarly to above, check the <code>/usr/lib/x86_64-linux-gnu/hdf5/serial/</code>
directory and look for file <code>libhdf5.a</code>. When found,
specify the location using -L flag or <code>CMakeLists.txt</code> by using
<syntaxhighlight lang="bash">
link_directories(/usr/lib/x86_64-linux-gnu/hdf5/serial/) # or your appropriate directory
</syntaxhighlight>
or fix the problem permanently by soft-linking
<syntaxhighlight lang="bash">
sudo ln -s /usr/lib/x86_64-linux-gnu/hdf5/serial/* /usr/lib
</syntaxhighlight>

==Linear Algebra==
We use [http://eigen.tuxfamily.org/ Eigen] as our matrix library. See
[http://eigen.tuxfamily.org/dox-devel/group__QuickRefPage.html here] for use
reference and documentation. For a quick transition from Matlab see
[http://eigen.tuxfamily.org/dox/AsciiQuickReference.txt here].

== Using Intel Math Kernel Library (MKL) ==
Install [https://software.intel.com/en-us/mkl Intel MKL] and take not of installation directories.

Proper include and link directories need to be set to use the Intel MKL.
<syntaxhighlight lang="cmake">
include_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/include) # change these to your installation path
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/mkl/lib/intel64)
link_directories(SYSTEM /opt/intel/compilers_and_libraries/linux/lib/intel64)
</syntaxhighlight>

Eigen has great support for MKL. You can see more detailed instructions [https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html on their website].
To use MKL for math operations, define <code>EIGEN_USE_MKL_VML</code> when compiling. You must also link
the appropriate libraries and define <code>MKL_LP64</code> for use on a 64bit system.
With parallelism enabled, the configuration looks as follows

<syntaxhighlight lang="cmake">
target_compile_options(my_target PRIVATE "-fopenmp")
target_compile_definitions(my_target PUBLIC EIGEN_USE_MKL_VML MKL_LP64)
target_link_libraries(my_target medusa mkl_intel_lp64 mkl_intel_thread mkl_core pthread iomp5)
</syntaxhighlight>

If you have Intel Parallel Studio installed this also enables you to use the Pardiso paralle direct sparse solver through its [https://eigen.tuxfamily.org/dox/group__PardisoSupport__Module.html Eigen interface].

== Using Intel C/C++ Compiler ==

In order to use Intel's compiler when compiling Medusa, you have several standard optionas for <code>cmake</code>.
Make sure compilers and installed and in your <code>PATH</code> by running <code>which icc</code>, which
should return something like <code>/opt/intel/bin/icc</code>.

You can define the compiler when *first* calling cmake like so
<syntaxhighlight lang="bash">
cmake .. -DCMAKE_C_COMPILER=$(which icc) -DCMAKE_CXX_COMPILER=$(which icpc)
</syntaxhighlight>
If this is not your first call, remove the <code>build</code> directory and start anew.

You can also set the <code>CXX</code> and <code>CC</code> bash variables. Before calling
<code> cmake </code> for the first time you have to export the following:
<syntaxhighlight lang="bash">
export CXX="icpc"
export CC="icc"
</syntaxhighlight>

Relaxation of the nodal distribution

2022-06-16T12:50:33Z

AndrejKP: Summing variable from n to i so the names don't clash.

To construct stable and reliable shape functions the support domains need to be non-degenerated <ref name="LeeLocal">Lee CK, Liu X, Fan SC. Local muliquadric approximation for solving boundary value problems. Comput Mech. 2003;30:395-409.</ref>, i.e. the distances between support nodes have to be balanced. Naturally, this condition is fulfilled in regular nodal distributions, but when working with complex geometries, the nodes have to be positioned accordingly. There are different algorithms designed to optimally fill the domain with different shapes <ref name="LohnerGeneral">Löhner R, Oñate E. A general advancing front technique for filling space with arbitrary objects. Int J Numer Meth Eng. 2004;61:1977-91.</ref><ref name="LiuNode">Liu Y, Nie Y, Zhang W, Wang L. Node placement method by bubble simulation and its application. CMES-Comp Model Eng. 2010;55:89.</ref>, see [[Positioning of computational nodes]] for our algorithms. Here, an intrinsic feature of the MLSM is used to take care of that problem. The goal is to minimize the overall support domain degeneration in order to attain stable numerical solution. In other words, a global optimization problem with the overall deformation of the local support domains acting as the cost function is tackled. We seek the global minimum by a local iterative approach. In each iteration, the computational nodes are translated according to the local derivative of the potential
\begin{equation}
\delta \b{p}\left( \b{p} \right)=-\sigma_{k}\sum\limits_{i=1}^{n}{\nabla }V\left( \b{p}-\b{p}_i \right)
\end{equation}
where $V, n, \delta \b{p}, \b{p}_{i}$ and $\sigma_{k}$ stand for the potential, number of support nodes, offset of the node, position of i-th support node and relaxation parameter, respectively. After offsets in all nodes are computed, the nodes are repositioned as
$\b{p}\leftarrow \b{p}+\delta \b{p}\left( \b{p} \right)$.
Presented iterative process procedure begins with positioning of boundary nodes, which is considered as the definition of the domain, and then followed by the positioning of internal nodes.

The <code>BasicRelax</code> engine supports two call types:
* with supplied distribution function <code>relax(domain, func)</code>, where it tries to satisfy the user supplied nodal density function <code>func</code>. This can be achieved only when there is the total number of domain nodes the same as integral of density function over the domain. If there is too much nodes a volatile relax might occur. If there is not enough nodes the relax might become lazy. The best use of this mode is in combination with <code>GeneralFill</code> and/or <code>GeneralSurfaceFill</code>.
* without distribution, where nodes always move towards less populated area. The relax magnitude is simply determined from Annealing factor and distance to the closest node. A simple and stable approach, however, note that this relax always converges towards uniformly distributed nodes.

Example of filling and relaxing a 2D domain can be found in below code snippet and Figures. Note the difference between relax without supplied distribution (<xr id="fig:relax_uniform"/>) and with supplied distribution (<xr id="fig:relax_non-uniform"/>). These examples are extreme, relaxation is usually used to convert an almost desirable node distribution into a desirable node distribution.
More examples can be found in main repository under tests.
<syntaxhighlight lang="c++" line>
// Define domain shape.
double r = 0.25;
BallShape<Vec2d> bs({0.5, 0.5}, r);

// Fill domain disretization with h.
auto h = [](Vec2d p){
return (0.005 + (p[0] - 0.5) * (p[0] - 0.5) / 2 + (p[1] - 0.5) * (p[1] - 0.5) / 2);
};
DomainDiscretization<Vec2d> domain = bs.discretizeBoundaryWithDensity(h);
GeneralFill<Vec2d> gf;
gf(domain, h);

// Relax domain to uniform distribution.
BasicRelax relax;
relax.iterations(1000).initialHeat(1).finalHeat(0).boundaryProjectionThreshold(0.75)
.projectionType(BasicRelax::PROJECT_IN_DIRECTION).numNeighbours(5);
relax(domain);
</syntaxhighlight>

<figure id="fig:relax_uniform">
[[File:relax_uniform.png|thumb|center|800px|<caption> Relaxation of a non-uniform node distribution to a uniform node distribution. </caption>]]
</figure>

<figure id="fig:relax_non-uniform">
[[File:relax_non-uniform.png|thumb|center|800px|<caption> Relaxation of a uniform node distribution to a non-uniform node distribution. </caption>]]
</figure>

=References=
<references/>

Radial basis function-generated finite differences (RBF-FD)

2022-06-16T12:49:06Z

AndrejKP: Added clarification to what PHS is since it wasn't mentioned until now.

This page describes the computation of RBF-FD weight augmented with polynomials. See also [[Computation of shape functions]] and [[Meshless Local Strong Form Method (MLSM)]] for a more general discussion.

== RBF-FD ==
$
\newcommand{\R}{\mathbb{R}}
\newcommand{\T}{\mathsf{T}}
\renewcommand{\L}{\mathcal{L}}
\renewcommand{\b}{\boldsymbol}
\newcommand{\n}{\b{n}}
\newcommand{\x}{\b{x}}
\newcommand{\w}{\b{w}}
\newcommand{\eps}{\varepsilon}
\newcommand{\lap}{\nabla^2}
\newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
$

Approximations of partial differential operators are the core
of strong form meshless procedures. Consider a partial differential operator $\L$
at a point $\x_c$. Approximation of $\L$ at a point $\x_c$
is sought using an ansatz
\begin{equation}
\label{eq:approx}
(\L u)(\x_c) \approx \sum_{i=1}^{n} w_i u(\x_i).
\end{equation}
Here $\x_i$ are the neighboring nodes of $\x_c$ which
constitute its stencil, $w_i$ are called stencil weights,
$n$ is the stencil size and $u$ is an arbitrary function.

This form of approximation is desirable, since operator $\L$ at point $\x_c$
is approximated by a linear functional $\w_\L (\x_c)$, assembled of weights $\w_i$
\begin{equation}
\label{eq:approx-vec}
\L|_{\x_c} \approx \w_\L (\x_c)^\T
\end{equation}
and the approximation is obtained using just a dot product with the function values
in neighboring nodes. The dependence of $\w_\L (\x_c)$ on $\L$ and $\x_c$ is often omitted
and written simply as $\w$.

To determine the unknown weights $\w$, equality of \eqref{eq:approx}
is enforced for a given set of functions. A natural choice are monomials, which
are also used in FDM, resulting in the Finite Point Method.

In the RBF-FD discretization the equality is satisfied for radial basis functions $\phi_j$.
Each $\phi_j$, for $j = 1, \ldots, n$ corresponds to one linear equation
\begin{equation}
\sum_{i=1}^{n} w_i \phi_j (\x_i) = (\L \phi_j)(\x_c)
\end{equation}
for unknowns $w_i$. Assembling these $n$ equations for into matrix form,
we obtain the following linear system:
\begin{equation} \label{eq:rbf-system}
\begin{bmatrix}
\phi(\|\x_1 - \x_1\|) &\cdots & \phi(\|\x_n - \x_1\|) \\
\vdots & \ddots & \vdots \\
\phi(\|\x_1 - \x_n\|) &\cdots & \phi(\|\x_n - \x_n\|)
\end{bmatrix}
\begin{bmatrix}
w_1 \\ \vdots \\ w_n
\end{bmatrix}
=
\begin{bmatrix}
(\L\phi(\|\x-\x_1\|))|_{\x=\x_c} \\
\vdots \\
(\L\phi(\|\x-\x_n\|))|_{\x=\x_c} \\
\end{bmatrix},
\end{equation}
where $\phi_j$ have been expanded for clarity.
The above system will be written more compactly as
\begin{equation} \label{eq:rbf-system-c}
A \w = \b \ell_\phi.
\end{equation}
The matrix $A$ is symmetric, and for some $\phi$ even positive
definite.<ref>Wendland, Holger. Scattered data approximation. Vol. 17. Cambridge university press, 2004.</ref> The [[Solving linear systems]] page has more details on numerical techniques for solving.

== Monomial augmentation ==
Using approximations that only contain RBF can lead to stability issues with
conditioning under refinement or failure to converge due to stagnation errors.
To improve accuracy and convergence, the approximation can be augmented with polynomials.

Let $p_1, \ldots, p_s$ be polynomials forming the basis of the space of $d$-dimensional
multivariate polynomials up to and including total degree $m$, with $s = \binom{m+d}{d}$.

Since enforcing exactness of \eqref{eq:approx}
for additional function would result in an overdetermined system, these additional constraints are
enforced by extending \eqref{eq:rbf-system-c}
as
\begin{equation} \label{eq:rbf-system-aug}
\begin{bmatrix}
A & P \\
P^\T & 0
\end{bmatrix}
\begin{bmatrix}
\w \\ \b \lambda
\end{bmatrix}
=
\begin{bmatrix}
\b \ell_{\phi} \\ \b \ell_{p}
\end{bmatrix}\!\!,\
P = \begin{bmatrix}
p_1(\x_1) & \cdots & p_s(\x_1) \\
\vdots & \ddots & \vdots \\
p_1(\x_n) & \cdots & p_s(\x_n) \\
\end{bmatrix}\!\!,\
\b \ell_p = \begin{bmatrix}
(\L p_1)|_{\x=\x_c} \\
\vdots \\
(\L p_s)|_{\x=\x_c} \\
\end{bmatrix}
\end{equation}
where $P$ is a $n \times s$ matrix of polynomials evaluated at stencil nodes and
$\b \ell_p$ is the vector of values assembled by applying considered operator $\L$ to
the polynomials at $\x_c$.

Note that the equation $P^\T \w = \b \ell_p$ contains exactly exactness constraints for
$p_j$. However, the introduction of parameters $\lambda_j$
causes \eqref{eq:approx} to not be exact for $\phi_i$ anymore. In fact, is was
shown to be equivalent to the following constrained
minimisation problem
\begin{equation}
\min_{\w} \left(\frac{1}{2} \w^\T A \w - \w^\T \b \ell_{\phi}\right), \text{ subject to }
P^\T \b
w = \ell_{p}
\end{equation}
and parameters $\b \lambda$ can be interpreted as Lagrangian multipliers.

Weights obtained by solving \eqref{eq:rbf-system-aug}
are taken as approximations of $\L$ at $\x_c$, while values $\b \lambda$ are discarded.

== Example ==
The $\ell_\infty$ error of solving Poisson equation on a unit ball using RBF-FD with PHS (Polyharmonic spline) $\phi(r) = r^3$ with various orders of monomial augmentation is shown below.
<ref name="SlakRef">Slak, J., & Stojanovič B. & Kosec, G. (2019). High order RBF-FD approximations with application to a scattering
problem. SpliTech 2019.</ref>

Exactness constraints for monomials often produce high-order methods, since local error of function approximation is small (imagine that the Taylor expansion is exact up to a high order).

[[File:poisson_augmentation_convergence.png|600px]]

=References=
<references/>

Integrators for time stepping

2022-06-16T12:45:03Z

AndrejKP: Typo fixes

This page describes how to solve ordinary differential equations numerically with examples from our library.

== Introduction and notation ==

We are solving an initial value problem, given as

$
\begin{align*}
\dot{y}(t) &= f(t, y) \\
y(t_0) &= y_0
\end{align*}
$

where $y$ is the unknown (possibly vector) function, $t_0$ is the start time, $f$ is the derivative (the functions we wish to integrate) and $y_0$ is the initial value of $y$.
Numerically, we usually choose a time step $\Delta t$ and integrate the function up to a certain time $t_{\max}$. Times of subsequent time steps are denoted with $t_i$ and function values with $y_i$.

The simplest method is explicit Euler's method:
$y_{n+1} = y_{n} + \Delta t f(t, y_n)$

== Explicit (single step) methods ==

A family of single step methods are explicit Runge-Kutta methods, which use intermediate derivative values to give a better estimation of value $y_{n+1}$ on the next step.

It is given by

$y_{n+1} = y_n + h \displaystyle \sum_{i=1}^s \beta_i k_i$

where

$
\begin{align*}
k_1 & = f(t_n, y_n), \\
k_2 & = f(t_n+\gamma_2h, y_n+h(\alpha_{21}k_1)), \\
k_3 & = f(t_n+\gamma_3h, y_n+h(\alpha_{31}k_1+\alpha_{32}k_2)), \\
& \ \ \vdots \\
k_s & = f(t_n+\gamma_sh, y_n+h(\alpha_{s1}k_1+\alpha_{s2}k_2+\cdots+\alpha_{s,s-1}k_{s-1})).
\end{align*}
$

To specify a particular method, one needs to provide the integer $s$ (the number of stages), and the coefficients $\alpha_{ij}$, $\beta_i$ and $\gamma_i$. This structure is known as the Butcher's tableau of the method:
$
\begin{array}{l|l}
\gamma & \alpha \\ \hline
& \beta^\T
\end{array}
$

First order method is the Euler's method above, while methods of very high order are available.
The most famous is RK4, the Runge Kutta method of fourth order.

$\begin{align*}
& {{k}_{1}}=f({{t}_{i}},{{y}_{i}}) \\
& {{k}_{2}}=f({{t}_{i}}+\Delta t/2,{{y}_{i}}+\Delta t{{k}_{1}}/2) \\
& {{k}_{3}}=f({{t}_{i}}+\Delta t/2,{{y}_{i}}+\Delta t{{k}_{2}}/2) \\
& {{k}_{4}}=f({{t}_{i}}+\Delta t,{{y}_{i}}+\Delta t{{k}_{3}}) \\
& {{y}_{i+1}}={{y}_{i}}+\Delta t/6({{k}_{1}}+2{{k}_{2}}+2{{k}_{3}}+{{k}_{4}}). \\
\end{align*}$

A more complete list can be found [https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods here].

Most common methods are implemented in the [http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/namespacemm_1_1integrators_1_1Explicit.html Explicit] namespace.

== Explicit multistep methods ==
These methods are called Adams–Bashforth methods, and need previous $s$ are used to estimate next value with a higher order. These method hope to gain computational efficiency by decreasing the number of function computations needed by using already computed values.

The explicit version is stated as $y_{n+s} = y_{n+s-1} + \sum_{i=1}^s b_i y_{n+s-i}$. The values $b_i$ are called weights.

The drawbacks of these mehtods are that another method of same order is needed to compute the initial time steps and that they can not be easily modified to work with adaptive time steps.

Most common methods are implemented in [http://www-e6.ijs.si/ParallelAndDistributedSystems/MeshlessMachine/technical_docs/html/namespacemm_1_1integrators_1_1ExplicitMultistep.html ExplicitMultistep] namespace.

== Adaptive methods ==
Adaptive methods estimate the error on each time step and descrease or increase the step as necessary.
Usually the have some addition parameters, such as $\Delta t_{\max}$ and $\Delta t_{\min}$ representing maximal and minimal allowed time steps and $\varepsilon$, represeting tolrance.

A family of adaptive methods can be derived from Runge Kutta methods above. Using a method of order $k$ and $k+1$, that differ only in weights, called $\beta$ and $\beta^\ast$, error can be estimated using [https://en.wikipedia.org/wiki/Richardson_extrapolation Richardson Extrapolation]:

$e_{n+1} = y_{n+1} - y^*_{n+1} = h\sum_{i=1}^s (\beta_i - \beta^*_i) k_i$.

If $\|e_{n+1}\| < \varepsilon$, time step can be increased, otherwise it is decreased. This can be achieved by scaling $\Delta t$ in accorance with $\varepsilon / \|e_{n+1}\|$.

These methods are implemented in Matlab as ode45 and similar. They can be found [https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Embedded_methods here], with a sample Matlab implementation of Cash-Karp method [https://github.com/jureslak/numerika-fmf/blob/master/ninde/hw2/CashKarp.m here].

These methods are not implemented yet, but will probably be added in the future.

== Usage ==
Ordinary ODE:
Solve $\dot{x} = y, \dot{y} = -x$ with $x(0) = 1, y(0) = 0$.

<syntaxhighlight lang="c++">
// define the function
std::function<Eigen::VectorXd(double, const Eigen::VectorXd&)> func =
[](double, const Eigen::VectorXd& y) {
Eigen::VectorXd r(2);
r(0) = -y(1);
r(1) = y(0);
return r;
};

double tmax = 2*M_PI;
double step = 0.1;
Eigen::VectorXd y0(2); y0 << 1.0, 0.0;
auto integrator = integrators::Explicit::RK4().solve(func, 0.0, tmax, step, y0); // clss representing the method and storing the function and initial conditions.

auto stepper = integrator_rk4.begin(); // iterator over the solver, integrating $f$ step by step

while (stepper_rk4) { // cast to bool tells us if we are done
++stepper_rk4;
// do something with step
}

// multiple steppers can exists at the same time
auto stepper2 = integrator.begin();

// we can also use range based for loop to integrate by time
for (auto& step : integrator) {
// do something with step
}
</syntaxhighlight>

Heat equation:

<syntaxhighlight lang="c++">
// after domain and operators are created
auto dv_dt = [&](double, const VecXd& y) {
Eigen::VectorXd der(y.size());
for (int c : interior) {
der[c] = op.lap(y, c);
}
for (int c : boundary) {
der[c] = 0;
}
return der;
};

auto integrator = integrators::Explicit::RK4().solve(dv_dt, 0.0, time, dt, T1); // choose your own integrator
for (auto& step : integrator) {
// do something
}
</syntaxhighlight>

[https://en.wikipedia.org/wiki/Lorenz_system Lorenz system:]

<syntaxhighlight lang="c++">

// System parameters
double sigma = 10.;
double rho = 28.;
double beta = 8./3.;

auto dy_dt = [&](double, const Vec3d& y) {
Vec3d der(3);
der(0) = sigma*(y(1)-y(0));
der(1) = y(0)*(rho-y(2))-y(1);
der(2) = y(0)*y(1)-beta*y(2);
return der;
};

double dt = 0.001; // Timestep
Vec3d y0 = {1, 1, 1}; // Initial condition

auto integrator = integrators::Explicit::RK4().solve(dy_dt, 0, 100, dt, y0);
for (auto& step : integrator) {
auto t = step.time();
auto y = step.value();
// output values
}
</syntaxhighlight>

<figure id="fig:lorenz_system">
[[File:lorenz.png|500px|thumb|upright=2|center|<caption>A sample solution of the Lorenz system for the values $\sigma=10$, $\rho=28$ and $\beta = 8/3$. The two trajectories (one in orange, the other in blue) correspond to two initial points that differ by a small amount in the $z$-coordinate.</caption>]]
</figure>

Ghost nodes (theory)

2022-06-16T12:40:37Z

AndrejKP: Typo fixes

<figure id="fig:tests">
[[File:ghost.png|417px|thumb|upright=2|alt=Ghost nodes schema.|<caption>Ghost node configuration at the domain boundary.</caption>]]
</figure>

Ghost nodes are a technique for discretizing (mostly) Neumann boundary conditions in PDEs. See the [[Ghost nodes|example.]]

== Introduction ==
Ghost nodes are a technique used for discretizing Neumann boundary conditions in FDM. To be able to use the central difference for first derivative, additional point, called ghost point, is introduced outside the domain boundary.
The unknown function value at the ghost node is added as a variable. At the boundary node, the Neumann condition is enforced, as well as the equation itself (two equations for two unknowns, the ghost and the boundary function value).

== Meshless setting (implicit) ==
Consider the problem $\mathcal{L}u = f$ with Neumann boundary conditions $\frac{\partial u}{\partial n} = g$ on some portion $\Gamma_1$ of the boundary and
Dirichlet conditions $u = u_0$ on portion $\Gamma_2$. Denote the number of internal nodes with $N_i$, number of Neumann boundary nodes with $N_n$ and number of Dirichlet boundary nodes with $N_d$.
The total number of nodes $N$ is equal to $N = N_i + N_n + N_d$ and so is the number of unknowns, representing solution values at these nodes.

For each Neumann node $p$, additional ''ghost'' or ''fictious'' nodes are placed outside the domain, as seen in the figure on the right. This is usually done
so that one ghost node is added for each boundary node, spaced $h$ away in the normal direction, where $h$ is the local nodal spacing.
This increases the number of nodes and unknowns by $N_n$.
Stencils and stencil weights are computed as before, and ghost nodes are included in the stencils. For each node $p_j$ denote the indices of its stencil nodes by $I_j$.
We will denote the computed stencil weights for operators with $w_\mathcal{L,j}$ and $w_{\frac{\partial }{\partial n},j}$.

We have the same equations as before:

* $N_i$ for interior nodes: $w_{\mathcal{L},j} \cdot \boldsymbol{u}_{I_j} = f_j$
* $N_d$ for Dirichlet nodes: $u_j = u_{0,j}$
* $N_n$ for Neumann nodes: $w_{\frac{\partial }{\partial n},j} \cdot \boldsymbol{u}_{I_j} = g_j$

Additional $N_n$ equations for ghost nodes are required. They are obtained by stating that the PDE itself must also hold in the boundary node, i.e., we add

* additional $N_n$ equations for Neumann nodes: $w_{\mathcal{L},j} \cdot \boldsymbol{u}_{I_j} = f_j$.

Note that no requirements about the PDE or BCs are made in the ghost nodes, they are simply involved as stencil nodes in the two equations for Neumann nodes (and possibly others).

All these equations are assembled in a sparse matrix which is solved and values obtained for ghost nodes can be neglected.

Additionally, more than one layer of ghost nodes can be added and additional equations on the boundary are needed. Sometimes, equations $\mathcal{L}^2 u = \mathcal{L}f$, etc... are used.

For an example of ghost nodes generation and use in Medusa see the [[Ghost nodes]] example.

=== Explicit version ===

When explicit time stepping is used, Neumann boundary conditions and ghost nodes are more difficult to implement.
We assume that the values of the interior in the $n$-th iteration are already known and denoted by $u_j^n$, as well as the
boundary and ghost values from the previous iteration $u_j^{n-1}$.

Consider a boundary node $p_j$ with an associated ghost node $p_{\gamma_j}$ and stencil indices $I_j = \{j, \gamma_j, i_1, \ldots, i_{n-2}\}$.
Denote the stencil weights at $p_j$ for $\mathcal{L}$ as $w_{\mathcal{L},j}$ and for $\frac{\partial }{\partial n}$ as $w_{\frac{\partial }{\partial n},j}$.

We have two equations at this node:

$$\sum_{i \in I_j} w_{\mathcal{L},j,i} u_i^{n} = f_j$$
$$\sum_{i \in I_j} w_{\frac{\partial }{\partial n},j, i} u_i^{n} = g_j$$

A system of these $2N_n$ equations can be solver in each iteration to obtain the next values of ghost and boundary nodes.

Furthermore, this system can be made explicit using single step of Gauss-Seidel iteration on a single step. We use the second equation to compute the
new value of $u_{\gamma_j}$ and the first to compute the new value of $u_j$.

For each ghost node $p_{\gamma_j}$:
$$u_{\gamma_j}^n = \frac{1}{w_{\frac{\partial }{\partial n},j,\gamma_j}} \left[ g_j - w_{\frac{\partial }{\partial n},j,j} u_j^{n-1} - \sum_{k=1}^{n-2}w_{\frac{\partial }{\partial n},j,i_k} u_{i_k}^{n|n-1} \right]$$

For each boundary node $p_j$:
$$u_{j}^n = \frac{1}{w_{\mathcal{L},j,\gamma_j}} \left[ f_j - w_{\mathcal{L},j,\gamma_j} u_{\gamma_j}^{n} - \sum_{k=1}^{n-2}w_{\mathcal{L},j,i_k} u_{i_k}^{n|n-1} \right]$$

The time index $n|n-1$ refers to either $n$ it the value of node $i_k$ has already been computed in this iteration (e.g. it is an internal node or already computed ghost or boundary node) or $n-1$ if the value is now known yet.

We are assuming that the stencil weight of the ghost node is nonzero in the discretization of the Neumann BC and that the weight of the boundary noe is nonzero in the discretization of $\mathcal{L}$. If this is not the case, the equations need to be rearranged differently.

K-d tree

2022-06-16T12:30:48Z

AndrejKP: Grammar

{{DISPLAYTITLE:''k''-d tree}}
In structured meshes, neighborhood relations are implicitly determined by the mapping from the physical to the logical space.
In unstructured mesh based approaches, support domains can be determined from the mesh connectivity relations.
However, in meshless methods, the stencil/support nodes are usually chosen as the nearest $n$ visible nodes.

In order to compute nearest neighbors effectively, fast spatial search structures have been developed, such as [https://en.wikipedia.org/wiki/Octree octree], [https://en.wikipedia.org/wiki/Ball_tree ball tree], [https://en.wikipedia.org/wiki/Cover_tree cover tree] and [https://en.wikipedia.org/wiki/K-d_tree $k$-d tree]. All of these structures first take a list of $N$ points and build a search structure, enabling subsequent queries for $n$ nearest neighbors to be answered in $O(n\log N)$ time.
The $k$-d tree structure is described below to give the reader a basic idea of how such spatial structures store data.

[[File:image_1avj3kcbe157a1ad610k81brmgsj9.png|800px|thumb|<caption> 2D tree example</caption>]]

A $k$-d 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 division of points is the 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$ coordinate 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.

Often, dynamic search structures are needed as well, which support node insertion and removal. In the beginning we were using [http://www.cs.umd.edu/~mount/ANN/ ANN] and [https://www.savarese.com/software/libssrckdtree/ libssrckdtree]
libraries for static and dynamic spatial search structures, respectively. Later, we tested four different publicly available C++ libraries that support at least static spatial search:

* [http://www.cs.umd.edu/~mount/ANN/ ANN]
* [https://www.savarese.com/software/libssrckdtree/ libssrckdtree]
* [http://spatial.sourceforge.net/ spatial]
* [https://github.com/jlblancoc/nanoflann/ nanoflann]

Two typical usage examples for these structures were tested. First, we tested a typical example of finding support domains of $n$ nodes for all nodes, which includes building a tree from a set of points and finding $n$ closest nodes for each point. The results are shown below.

[[File:kdtree_find_support_time.png]]

Second, we tested how the structures performed during the node generation process, as described in [[Positioning of computational nodes]] page. Results are shown in figure below.

[[File:kdtree_pds_time.png]]

Besides speed, which was the main factor, we also considered the number of dependencies a library has and how general it is regarding distances and data types.
Overall, we decided to use the [https://github.com/jlblancoc/nanoflann/ nanoflann] library in [http://e6.ijs.si/medusa/ medusa]. Our wrappers are available
as [http://e6.ijs.si/medusa/docs/html/classmm_1_1KDTree.html KDTree] and [http://e6.ijs.si/medusa/docs/html/classmm_1_1KDTreeMutable.html KDTreeMutable].

The results and reasoning is described more accurately in the technical report [[:File:report.pdf]].



== References ==
* Trobec R., Kosec G., Parallel scientific computing: theory, algorithms, and applications of mesh based and meshless methods: Springer; 2015.

Refinement of the nodal distribution

2022-06-16T12:29:04Z

AndrejKP: Typo fix. Their is -> there is

Here we consider possible meshless refinement algorithms (sometimes also called adaptive cloud refinement). The refinement mechanisms we have so far studied include:
* refinement based on closest node distance
* refinement based on averaged (inter-)node distance
* refinement based on half-links

Here we only want to compare the quality of the refined grids and have not tied the refinement algorithm with a error indicator, thus we only study the node insertion process by refining the whole grid.

The refinement routine takes a range of nodes (e.g. a subregion of the domain) together with the refinement parameters and generates new nodes around the old ones. Special care must be taken with refinement of the boundary nodes. Points have to be selected on the actual boundary either analytically considering the geometry or with a numerical root finder such as bisection.

==Problem description==

To compare the node refinement mechanisms we study the process of reaction-diffusion in an infinite cylindrical catalyst pellet (infinite in the $z$-dimension). Since the pellet is infinite in one dimension this problem simplifies to a 2D problem (in the $xy$-plane). For a catalyst pellet of radius $R$ centered at $(x,y) = (0,0)$ and the reactant undergoing a first order reaction we must solve the equation
\begin{equation}
\b{\nabla}^2 C - {M_T}^2 C = 0,
\end{equation}
where $C$ is the concentration of the reactant, $M_T = R\sqrt{k/D}$ is known as Thiele's modulus and $k$ and $D$ represent the reaction rate constant and diffusivity of the reacting species. The boundary conditions for this problem is \[C(R) = C_s.\] The analytical solution can be found easily using cylindrical coordinates and is given by
\begin{equation}
\frac{C(r)}{C_S} = \frac{I_0(r M_T)}{I_0(R M_T)},
\end{equation}
where $I_0(r)$ is the modified Bessel function of first kind (this function is available in the library Boost as well as scripting languages such as Python or MATLAB). The conversion from cartesian to cylindrical coordinates is given by \[r = \sqrt{x^2+y^2}.\]

==Error indicators==

To compare the quality of the refined meshes for the described problem case we look at different error criteria including the max norm $L_\infty$ defined as
\begin{equation}
L_\infty = \mathrm{max}_i \left|C_i^\mathrm{numerical} - C_i^\mathrm{analytical}\right|,
\end{equation}
the $L_2$ norm per node defined as
\begin{equation}
\bar{L_2} = \frac{\sqrt{\sum^N_{i = 1}\left(C_i^\mathrm{numerical} - C_i^\mathrm{analytical}\right)^2}}{N},
\end{equation}
where $N$ is the number of nodes (and pertinent equations) in the domain.

We also measure the number of iterations required by the sparse BiCGSTAB solver to reach convergence and the estimated error of solving the system of equations.

== Closest node ==

For a given node $\b{x}_0 = (x_0,y_0)$:

# find the closest node $\b{x}_1 = (x_1,y_1)$
# calculate the half distance between the two nodes \[d = |\b{x}_1 - \b{x}_0|/2\]
# randomly select up to 6 (The case of 6 nodes is the limit since it produces a regular hexagon. In practice this never occurs due to the "monte carlo" node selection procedure.) new nodes on the circle with center $\b{x}_0$ and radius $d$ and simultaneously make sure there is a minimal inter-nodal distance $d$ between the new nodes.

For boundary points we first select 2 points that intersect with the boundary of the domain and only then points lying inside the domain. Due to geometrical constraints boundary points will usually end up with 3 new nodes (in case of straight boundaries we could end up with 4, which would be the previously discussed hexagon limit).

<figure id="fig:node_refinement_1">
[[File:closest_node.png|1000px|thumb|center|<caption>Refinement based on closest node approach (initial unrefined grid is on the left). In the second refinement step an erroneous point has appeared from an internal point that was too close to the boundary. Also noticable is clustering of points on the boundary.</caption>]]
</figure>

== Average radius ==

Input parameters: $f$ and $l_s$

For a given node $\b{x}_0$:
# find the $l_s$ (an integer from 1 to 7) closest nodes
# calculate the average distance $\bar{d}$ to the $l_s$ closest nodes
# randomly select up to 6 new nodes on the circle with center $\b{x}_0$ and radius $f\cdot\bar{d}$ where $f$ is the radius fraction that lies between 0.2 (leads to clustering) and 0.8. Only allow nodes that are separated by the distance $f \cdot \bar{d}$.

''(note that in case $l_s = 1$ and $f = 0.5$ the average radius mechanism becomes equal to the closest node refinement approach described above)''

<figure id="fig:node_refinement_2">
[[File:new_average_radius.png|1000px|thumb|center|<caption>Refinement based on average radius approach (initial unrefined grid is on the left). The parameters are $l_s = 5$ closest nodes in average radius calculation and points placed at radius fraction $f = 0.5$. </caption>]]
</figure>

== Half-links ==

''Warning: this is still a prototype''

Input parameters: $l_s$, $d_m$

For a given node $\b{x}_0$:
# find the $l_s$ (an integer from 1 to 7) closest nodes $\b{x}_i$
# select new nodes in the middle of the segments $\b{x}_i - \b{x}_0$ only allowing points that are separated by the minimal distance $d_m$

''(note also that in the 1D case the half-link and closest radius approach become the same)''

The minimal distance $d_m$ is chosen as a fraction of the distance to the closest link, e.g. $d_m = f d$, where $f$ is the provided fraction and $d$ is the distance to the closest link.

<figure id="fig:node_refinement_3">
[[File:half_link.png|1000px|thumb|center|<caption>Refinement based on half links (initial unrefined grid is on the left). The parameters are $l_s = 6$ and $d_m = 0.4 d$, where $d$ is the distance to the closest link.</caption>]]
</figure>

<figure id="fig:node_refinement_3a">
[[File:refineWithRelax.png|1000px|thumb|center|<caption>Refinement based on half links with additional 10 step relax after refinement</caption>]]
</figure>

After experimentation we noticed there are some inconsistencies when trying to refine structured point sets with this approach. The reason for these inconsistencies is that the boundary and internal points have a different number of "natural neighbours". For example in 2D on a square grid, the internal points have 8 neighbours, while boundary points have 5 neighbours. If we choose higher numbers e.g 9 links for an internal node, the 9th node might be any of the 4 nodes one shell further out that only differ at machine precision.

The figures below show some preliminary results of refinement based on half-links. For the circle domain relaxation was applied after the refinement.

<figure id="fig:square_refinement">
[[File:dc_field.png|1000px|thumb|center|<caption> Thermal diffusion in (convective) flow at a stagnation point (bottom left corner).</caption>]]
</figure>

After the second refinement of the corner, the solver had difficulty converging to the solution. This was the result of a fixed size shape parameter in the shape functions of the node points. The shape functions have to be tailored to the local characteristic distance in the point set.

<figure id="fig:circle_refinement">
[[File:ref_circle.png|1000px|thumb|center|<caption>Reaction-diffusion in a cylinder catalyst. Two successive refinements have been applied for $r > 0.5$ and $r > 0.8$, where $r$ is the radial coordinate. The cylinder radius $R = 1$.$</caption>]]
</figure>

== Hybrid approach ==

A hybrid might give better distributions of the refined points. The half-link approach performs well at the boundaries while the distance approach gives less regular internal distributions. In any case it is suggested to perform a few more relaxation steps to equilibrate the mesh.

<figure id="fig:node_refinement_4">
[[File:hybrid_refine.png|1000px|thumb|center|<caption>Refinement based on half links at the boundaries and closest distances for the internal nodes (initial unrefined grid is on the left). for the boundary nodes the parameters are $l_s = 7$ and $d_m = 0.5 d$, where $d$ is the distance to the closest link.</caption>]]
</figure>