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

Realistic 3D models

2020-11-19T10:31:31Z

Jureslak: /* Working with 3D models in Medusa */

Medusa is not meant to deal with complex 3D geometry - we support balls, cuboids, and their differences, unions, translations and rotations. More general 2D shapes can be described by polygons, which are also included in Medusa. For more general 3D objects, the geometric algorithms become more complex and it makes sense to leverage third-party libraries.

== Formats, software, and other resources ==

A good library in c++ is CGAL ([https://www.cgal.org/]), that we can use to work with surface meshes: [https://doc.cgal.org/latest/Surface_mesh/index.html https://doc.cgal.org/latest/Surface_mesh/index.html]
CGAL can usually be installed with <code>apt-get install libcgal-dev</code> or similar. The Qt dependency is not necessary, but can be useful for easier visualization without exporting to an external program.

Polyhedrons are solid bodies with polygonal faces, and their surface forms a closed mesh. CGAL can work also with surfaces meshes that are not closed, but for PDE solving in 3D bodies, we assume them to be closed (or watertight) and triangular. Polyhedron data is usually stored in a separate file with many different formats available. Some of the formats actually store the polyhedron data (including how faces are connected) and some only store the faces. By default, CGAL reads OFF files [https://en.wikipedia.org/wiki/OFF_(file_format)], and these files can also be easily viewd by <code>geomview</code> [http://www.geomview.org/], also trivially installable on most Linux systems. OFF files include an option to specify the color of each face, which can be used to identify different surface regions.

Other formats include STL, PLY and OBJ, and other proprietary formats. Some of these formats are more appropriate than others. Notably STL, while very common, can be problematic, as it does not include connections between faces, bit only a list of triangles, from which the polyhedron must first be reconstructed. This can be done directly in CGAL by using <code>CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh</code>, through an intermediate program, like Meshlab https://www.meshlab.net/, or through an online converter, such as https://cadcook.com/ply-off/. Note that it is not necessary that any of these conversions work, since STL file can represent an object with holes in the surface mesh or other degenerated objects.
PLY files on the other hand do represent surface meshes, and can be read in CGAL (or easily converted to OFF) [https://cadcook.com/ply-off]. Similar hold for OBJ files, although, these files can contain many other types of objects besides polygon meshes. CGAL includes readers and writers for OFF, PLY, OBJ and STL formats [https://github.com/CGAL/cgal/tree/master/Polyhedron_IO/include/CGAL/IO], and the default <code>cin >> </code> operator expects an OFF file. Further demos for polyhedron IO are available in the CGAL repo [https://github.com/CGAL/cgal/tree/master/Polyhedron_IO].

PLY, OFF and STL files can also be read with Matlab, but require (short) external functions, available at [https://www.mathworks.com/matlabcentral/fileexchange/46540-geometry-processing-package], see Function/IO. These return a list of vertices and faces, that can be plotted with <code>patch</code> or <code>trisurf</code> after constructing the triangulation.

Obtaining models can be somewhat difficult, as this files are often used to store general 3D objects for rendering. Here is a list of common 3D models (not all of which are polyhedrons): [https://github.com/alecjacobson/common-3d-test-models https://github.com/alecjacobson/common-3d-test-models]. Specifically, an altered Stanford bunny model (made watertight) in OFF format can be found here [https://github.com/OpenGP/OpenGP/blob/master/data/bunny.off]. Googling with keywords "polyhedron", "surface mesh", "watertight", "closed", "OFF" or "PLY" might help.

General models (often STL files) can be downloaded from [https://grabcad.com/ https://grabcad.com/] (after registering). Usually it is simplest to inspect them with Meshlab and also use it to save files as OFF.

This shows the <code>bunny.off</code> file referenced above open in Meshlab.

[[File:bunny_meshlab.png|400px]]

Meshlab can also be used to fill holes, clean meshes and color faces. See the Filters / Remeshing or the Filters / Color menu. Meshlab also colors vertices, and when saving, one can choose to save only the face colors. Saving also the vetor files with cause the file to be saved as COFF format, which also supports vertex colors.

Below is an example of colored bunny. The colors can be used to identify surface regions where different boundary conditions should be applied.

[[File:bunny-colored.png|400px]]

== Working with 3D models in Medusa ==

We assume that an OFF file has been successfully obtained, in this example we will use the <code>bunny.off</code> file from [https://github.com/OpenGP/OpenGP/blob/master/data/bunny.off https://github.com/OpenGP/OpenGP/blob/master/data/bunny.off] that was colored as the picture above. To work with Medusa, we can make a generic <code>PolyhedronShape</code>, that can also be used with other existing shapes (translations, boolean operations, etc).

This is done by inheriting from <code>DomainShape</code> and defining the required methods. Notably, CGAL is use to get a proper implementation of <code>contains</code> method, that is needed for filling the domain. CGAL is also used to get the bounding box, and to store the surface mesh in an easily traversable data structure. Note that the polyhedron, represented in the OFF file is considered to be the exact domain, which is then discretize with meshless techniques. The mesh is used only to define the domain, and is considered fixed and final. Discretizations can (and usually should) be smaller than.

Below is a picture of the bunny OFF model with surface nodes of various types, as defined by the face colors, along with the surface normals.

[[File:bunny-surface.png|400px]]

The <code>PolytopeShape</code> (a dimension independent version of <code>PolyhedronShape</code> and <code>PolygonShape</code>) can be used to work with OFF files. Otherwise, the solution procedure is as usual.

<syntaxhighlight lang="c++" line>
#include <medusa/bits/domains/PolyhedronShape.hpp> // This must be included separately.
...

int main() {
std::string filename = "bunny.off";
PolytopeShape<Vec3d> bunny = PolytopeShape<Vec3d>::fromOFF(filename);

double dx = 2.5;
DomainDiscretization<Vec3d> domain = bunny.discretizeBoundaryWithStep(dx);
GeneralFill<Vec3d> fill; fill.seed(0);
fill(domain, dx);

... // The usual.

for (int i : domain.interior()) {
0.1*op.grad(i, -1) + 2.0*op.lap(i) = -1.0;
}
for (int i : domain.boundary()) {
op.value(i) = 0.0;
}

Eigen::BiCGSTAB<decltype(M), Eigen::IncompleteLUT<double>> solver;
solver.preconditioner().setDroptol(1e-4);
solver.preconditioner().setFillfactor(20);
solver.compute(M);
ScalarFieldd u = solver.solve(rhs);

}
</syntaxhighlight>

The solution is shown below:

[[File:bunny-sol2.png|496px]]

2020-11-12T11:22:53Z

Jureslak: /* Working with 3D models in Medusa */

Medusa is not meant to deal with complex 3D geometry - we support balls, cuboids, and their differences, unions, translations and rotations. More general 2D shapes can be described by polygons, which are also included in Medusa. For more general 3D objects, the geometric algorithms become more complex and it makes sense to leverage third-party libraries.

== Formats, software, and other resources ==

A good library in c++ is CGAL ([https://www.cgal.org/]), that we can use to work with surface meshes: [https://doc.cgal.org/latest/Surface_mesh/index.html https://doc.cgal.org/latest/Surface_mesh/index.html]
CGAL can usually be installed with <code>apt-get install libcgal-dev</code> or similar. The Qt dependency is not necessary, but can be useful for easier visualization without exporting to an external program.

Polyhedrons are solid bodies with polygonal faces, and their surface forms a closed mesh. CGAL can work also with surfaces meshes that are not closed, but for PDE solving in 3D bodies, we assume them to be closed (or watertight) and triangular. Polyhedron data is usually stored in a separate file with many different formats available. Some of the formats actually store the polyhedron data (including how faces are connected) and some only store the faces. By default, CGAL reads OFF files [https://en.wikipedia.org/wiki/OFF_(file_format)], and these files can also be easily viewd by <code>geomview</code> [http://www.geomview.org/], also trivially installable on most Linux systems. OFF files include an option to specify the color of each face, which can be used to identify different surface regions.

Other formats include STL, PLY and OBJ, and other proprietary formats. Some of these formats are more appropriate than others. Notably STL, while very common, can be problematic, as it does not include connections between faces, bit only a list of triangles, from which the polyhedron must first be reconstructed. This can be done directly in CGAL by using <code>CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh</code>, through an intermediate program, like Meshlab https://www.meshlab.net/, or through an online converter, such as https://cadcook.com/ply-off/. Note that it is not necessary that any of these conversions work, since STL file can represent an object with holes in the surface mesh or other degenerated objects.
PLY files on the other hand do represent surface meshes, and can be read in CGAL (or easily converted to OFF) [https://cadcook.com/ply-off]. Similar hold for OBJ files, although, these files can contain many other types of objects besides polygon meshes. CGAL includes readers and writers for OFF, PLY, OBJ and STL formats [https://github.com/CGAL/cgal/tree/master/Polyhedron_IO/include/CGAL/IO], and the default <code>cin >> </code> operator expects an OFF file. Further demos for polyhedron IO are available in the CGAL repo [https://github.com/CGAL/cgal/tree/master/Polyhedron_IO].

PLY, OFF and STL files can also be read with Matlab, but require (short) external functions, available at [https://www.mathworks.com/matlabcentral/fileexchange/46540-geometry-processing-package], see Function/IO. These return a list of vertices and faces, that can be plotted with <code>patch</code> or <code>trisurf</code> after constructing the triangulation.

Obtaining models can be somewhat difficult, as this files are often used to store general 3D objects for rendering. Here is a list of common 3D models (not all of which are polyhedrons): [https://github.com/alecjacobson/common-3d-test-models https://github.com/alecjacobson/common-3d-test-models]. Specifically, an altered Stanford bunny model (made watertight) in OFF format can be found here [https://github.com/OpenGP/OpenGP/blob/master/data/bunny.off]. Googling with keywords "polyhedron", "surface mesh", "watertight", "closed", "OFF" or "PLY" might help.

General models (often STL files) can be downloaded from [https://grabcad.com/ https://grabcad.com/] (after registering). Usually it is simplest to inspect them with Meshlab and also use it to save files as OFF.

This shows the <code>bunny.off</code> file referenced above open in Meshlab.

[[File:bunny_meshlab.png|400px]]

Meshlab can also be used to fill holes, clean meshes and color faces.

== Working with 3D models in Medusa ==

We assume that an OFF file has been successfully obtained, in this example we will use the <code>bunny.off</code> file from [https://github.com/OpenGP/OpenGP/blob/master/data/bunny.off https://github.com/OpenGP/OpenGP/blob/master/data/bunny.off]. To work with Medusa, we can make a generic <code>PolyhedronShape</code>, that can also be used with other existing shapes (translations, boolean operations, etc).

This is done by inheriting from <code>DomainShape</code> and defining the required methods. Notably, CGAL is use to get a proper implementation of <code>contains</code> method, that is needed for filling the domain. CGAL is also used to get the bounding box, and to store the surface mesh in an easily traversable data structure.

<syntaxhighlight lang="c++" line>
#include <>

int main() {
}
</syntaxhighlight>

Realistic 3D models

2020-11-12T11:01:41Z

Jureslak:

Realistic 3D models

2020-11-12T11:01:24Z

Jureslak:

File:Bunny meshlab.png

2020-11-12T10:31:09Z

2020-11-09T13:51:25Z

Jureslak: /* Building blocks */

'''Welcome to the Medusa wiki. To visit the main website, go to [http://e6.ijs.si/medusa/ http://e6.ijs.si/medusa/].'''

In [http://e6.ijs.si/ParallelAndDistributedSystems/ Parallel and Distributed Systems Laboratory] we are working on a C++ library that is first and foremost focused on tools for solving Partial Differential Equations by meshless methods. The basic idea is to create generic codes for tools that are needed for solving not only PDEs but many other problems, e.g. Moving Least Squares approximation, $k$-d tree, domain generation engines, etc.
We call this open source meshless project [http://e6.ijs.si/medusa/ Medusa: Coordinate Free Meshless Method implementation (MM)].

Technical details about code and examples can be found on our [http://e6.ijs.si/medusa/docs/ documentation page] and [https://gitlab.com/e62Lab/medusa Gitlab repository]. [[File:C.png|100px||link=https://gitlab.com/e62Lab/medusa|alt=Alt text|code]] [[File:doxygen.png|100px|link=http://e6.ijs.si/medusa/docs/|alt=Alt text|Documentation page]]

This wiki site is meant for more relaxed discussions about general principles, possible and already implemented applications, preliminary analyses, etc.
Note, that there are many grammatical mistakes, typos, stupid sentences, etc. This wiki is meant for quick information exchange and therefore we do not invest a lot of energy into styling :).

== Documentation ==
* [https://gitlab.com/e62Lab/medusa Code on Gitlab]
* [[How to build | Installation and building]]
* [[Including this library in your project | Including this library in your project]]
* [[Testing | Running tests]]
* [http://e6.ijs.si/medusa/docs/ Technical documentation]
* [[Coding style | Coding style]]
* [[Wiki editing guide | Wiki editing and backup guide]]

== Examples ==
In this section we present exact examples. Each of the below solutions can be found also in in the repository under examples. More explanation about the physical background and solution procedure can be found in following sections.
* [[Philosophy of examples and how to run them]]
* [[Poisson's equation]]
* [[Heat equation]]
* [[Linear elasticity]]
* [[Complex-valued problems]]
* [[Coupled domains]]
* [[Parametric domains | Parametric domains – Curved surface with variable density]]
* [[NURBS domains | Domains modeled with non-uniform rational basis spline's (NURBS)]]
* [[Customization]]
* [[Ghost nodes]]
* [[Electromagnetic scattering]]
* [[Schrödinger equation]]
* [[Wave equation]]
* [[Cahn-Hilliard equation]]
* [[Non-Newtonian_fluid_example | Non-Newtonian fluid]]
* [[Meshless Lattice Boltzmann method]]

== Building blocks ==
Medusa is modular coordinate-free parallel implementation of a numerical framework designed, but not limited to, for solving PDEs. In this section we present main modules of the library that can be also used as a standalone tools.
* [[Positioning of computational nodes]]
* [[Relaxation of the nodal distribution]]
* [[Refinement of the nodal distribution]]
* [[k-d tree|''k''-d tree]] and other spatial search structures
* Solving [[Solving linear systems | linear systems]], [[Solving overdetermined systems | overdetermined]] and [[Solving underdetermined systems | underdetermined]]
* [[Weighted Least Squares (WLS)]]
* [[Computation of shape functions]]
* [[Meshless Local Strong Form Method (MLSM)]]
* [[Radial basis function-generated finite differences (RBF-FD)]]
* [[Ghost nodes (theory)]]
* [[Integrators for time stepping]]
* [[RBF Interpolation]]
* [[Complex 3D models]]

== Discussions / Applications ==
This section is meant for general discussion about the physical background of the examples, the solution procedures, various applications, etc. Note, that code snippets presented in discussion might not reflect the actual state of Medusa.
* Basic PDE solutions
** [[Convection Diffusion equation | Convection Diffusion equation]]
**[[Wave equation application]]
* [[Adaptivity]]
* [[Solid Mechanics]]
** [[Point contact]]
** [[Hertzian contact]]
** [[Cantilever beam]]
** [[Fretting fatigue case]]
* [[Fluid Mechanics]]
** [[Lid driven cavity]]
** [[de Vahl Davis natural convection test]]
** [[Natural convection in 3D irregular domain]]
** [[Natural convection from heated cylinder]]
** [[Natural convection between concentric cylinders]]
** [[Non-Newtonian fluid]]
* [[Computational electromagnetics]]
** [[Triple dielectric step in 1D]]
** [[Scattering from an infinite cylinder]]
** [[Point source near an anisotropic lens]]
* Other applications
** [[Attenuation due to liquid water content in the atmosphere|Attenuation of a satellite communication]]
** [[Heart rate variability detection]]
** [[Bioheat equation]]

== Performance analyses ==
* [[Execution on Intel® Xeon Phi™ co-processor]]
* [[1D MLSM and FDM comparison]]
* [[:File:tech_report.pdf|Execution overheads due to clumsy types::technical report]] [[File:pdf-file.gif]]
* [[Solving sparse systems]]
* [[Eigen paralelization]]

== Last changes ==
<news unique=1 limit = 5>
*{{{timeanddate}}} :: {{{title}}}

</news>

== Miscellaneous ==
* FAQ - [[Frequently asked questions]].
* [[List of wiki contributors]]
* List of library contributors: [http://e6.ijs.si/medusa/about#about-contributors See the official website]

== References ==
* Slak J., Kosec G. Adaptive radial basis function-generated finite differences method for contact problems. International journal for numerical methods in engineering, ISSN 0029-5981 [http://www-e6.ijs.si/ParallelAndDistributedSystems/pdf/32230439.pdf manuscript]
* Slak J., Kosec G.; Refined meshless local strong form solution of Cauchy-Navier equation on an irregular domain. Engineering analysis with boundary elements. 2018;11 ; [http://comms.ijs.si/~gkosec/data/papers/31107623.pdf manuscript]
* Depolli, M., Kosec, G., Assessment of differential evolution for multi-objective optimization in a natural convection problem solved by a local meshless method. Engineering optimization, 2017, vol. 49, no. 4, pp. 675-692 ;[http://comms.ijs.si/~gkosec/data/papers/29639719.pdf manuscript]
* Kosec G., A local numerical solution of a fluid-flow problem on an irregular domain. Advances in engineering software. 2016;7 ; [29512743] ; [http://comms.ijs.si/~gkosec/data/papers/29512743.pdf manuscript]
* Kosec G., Trobec R., Simulation of semiconductor devices with a local numerical approach. Engineering analysis with boundary elements. 2015;69-75; [27912487] ; [http://comms.ijs.si/~gkosec/data/papers/27912487.pdf manuscript]
* Kosec G., Šarler B., Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method. Engineering analysis with boundary elements. 2014;36-44; [http://comms.ijs.si/~gkosec/data/papers/3218939.pdf manuscript]
* Kosec G., Depolli M., Rashkovska A., Trobec R., Super linear speedup in a local parallel meshless solution of thermo-fluid problem. Computers & Structures. 2014;133:30-38; [http://comms.ijs.si/~gkosec/data/papers/27339815.pdf manuscript]
* Kosec G., Zinterhof P., Local strong form meshless method on multiple Graphics Processing Units. Computer modeling in engineering & sciences. 2013;91:377-396; [http://comms.ijs.si/~gkosec/data/papers/26785063.pdf manuscript]
* Kosec G., Šarler B., H-adaptive local radial basis function collocation meshless method. Computers, materials & continua. 2011;26:227-253; [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerBurgers.pdf manuscript]
* Trobec R., Kosec G., Šterk M., Šarler B., Comparison of local weak and strong form meshless methods for 2-D diffusion equation. Engineering analysis with boundary elements. 2012;36:310-321; [http://comms.ijs.si/~gkosec/data/papers/EABE2499.pdf manuscript]
* Kosec G, Zaloznik M, Sarler B, Combeau H. A Meshless Approach Towards Solution of Macrosegregation Phenomena. CMC: Computers, Materials, & Continua. 2011;580:1-27 [http://comms.ijs.si/~gkosec/data/papers/KosecZaloznikSarlerCombeauSegregation.pdf manuscript]
* Kosec G, Sarler B. Solution of thermo-fluid problems by collocation with local pressure correction. International Journal of Numerical Methods for Heat & Fluid Flow. 2008;18:868-82 [http://comms.ijs.si/~gkosec/data/papers/KosecSarlerNS2008.pdf manuscript]
* Trobec R., Kosec G., Parallel Scientific Computing, ISBN: 978-3-319-17072-5 (Print) 978-3-319-17073-2.
* Slak, J., Kosec, G.. Detection of heart rate variability from a wearable differential ECG device., MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938, pp 450-455.
* Kolman, M., Kosec, G. Correlation between attenuation of 20 GHz satellite communication link and liquid water content in the atmosphere. MIPRO 2016, 39th International Convention, 2016, Opatija, Croatia, ISSN 1847-3938. pp. 308-313.

==Related pages==
* http://e6.ijs.si/ParallelAndDistributedSystems/#!NumericalMethods
* http://e6.ijs.si/ParallelAndDistributedSystems/#!utils
* http://e6.ijs.si/ParallelAndDistributedSystems/#!NUMA

Computation of shape functions

2020-08-29T09:26:34Z

Jureslak: /* Another view on derivation of shape functions */

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 parial 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
\[
\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 depend 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.

Computation of shape functions

2020-08-29T09:24:41Z

Jureslak: /* Another view on derivation of shape functions */

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 parial 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
\[
\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 depend 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) (\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.

How to build

2020-07-06T19:09:54Z

Jureslak: /* Linker errors */

=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
</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!

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
</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>

How to build

2020-07-06T19:09:09Z

Jureslak:

=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
</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!

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 (without <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
</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>

How to build

2020-07-06T18:56:11Z

Jureslak:

=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
</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!

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 ==

== 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
</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>