.. nlfem documentation master file, created by
   sphinx-quickstart on Thu May  6 17:59:33 2021.
   You can adapt this file completely to your liking, but it should at least
   contain the root `toctree` directive.

Welcome to nlfem's documentation!
=================================
.. toctree::
   :maxdepth: 2
   :caption: Contents:

.. include:: ../../README.rst

-----------------------------------------------------------------

Quick Start
===========

To test the rates for the constant kernel run

::

  cd nlfem/examples/test2D
  python3 computeRates2D.py -s 4 -f testConfConstant

Run a test with other configurations by indicating the file e.g. ``-f testConfFull``. The option ``-s`` indicates the number of refinement steps.

Usage
======


Let :math:`\Omega \subset \mathbb{R}^{d}, (d=1,2,3)` be a bounded domain and :math:`\Gamma^D \subset \Omega^c`, where we allow that :math:`\Gamma^D = \emptyset`.
The package nlfem allows to assemble stiffness matrices which discretize the nonlocal operator

.. math::

    -\mathcal{L}_\delta \mathbf{u}(\mathbf{x}) = 2\int_{\Omega \cup \Gamma^D}
	 \mathbb{1}_{B_{\delta}(\mathbf{x})}(\mathbf{y})
	\left(\mathbf{C}(\mathbf{x}, \mathbf{y})\mathbf{u}(\mathbf{x})
	- \mathbf{C}(\mathbf{y}, \mathbf{x})\mathbf{u}(\mathbf{y}) \right)
	
where 
:math:`\mathbf{C} \in \mathbb{R}^{n\times n}, (n \geq 1)` is  a possibly matrix-valued integral kernel and :math:`B_{\delta}(\mathbf{x})` a :math:`p`-norm ball 
(:math:`p \in [1, \infty]`) of radius :math:`\delta`. Details about the mathematical background of nlfem are presented in [Klar2022]_.


**Example** We exemplify the usage of nlfem by solving a scalar, nonlocal Dirichlet-type problem

.. math::

	\begin{cases}
		-\mathcal{L}_\delta u(\mathbf{x}) = f_{\Omega}(\mathbf{x}) & \mathbf{x} \in \Omega\\
		u(\mathbf{x}) = g(\mathbf{x})					   & \mathbf{x} \in \Gamma^D
	\end{cases}

where :math:`\Omega` by a 2d-disk of radius 0.9 with a nonlocal boundary :math:`\Gamma^D` of width :math:`\delta=0.1`.  

	
Finite element mesh
--------------------

The construction requires a finite element mesh in a 1d, 2d, or 3d domain. The mesh is characterized by its ``elements`` and ``vertices``. The ``elements`` are given by a ``numpy.ndarray`` of datatype ``np.int_`` and shape ``(nE, d+1)``, where ``nE`` is the number of elements and ``d`` the dimension of the domain. The ``vertices`` are a numpy array of floats with shape ``(nV, d)``, where ``nV`` is the number of vertices.

The domain described by the the above arrays needs to be divided into different parts according to their purpose. To that end we define a numpy ``nd.array`` called ``elementLabels`` of type ``np.int_`` and length ``nE``. It assigns a label to each element. Negative elements indicate the Dirichlet boundary, any positive element is considered a part of the domain and zero labeled elements are not part of the integration. 

**Example** The finite element mesh 
(that is  ``vertices`` and ``elements``) of the 
required form can be obtained for example from `gmsh <http://gmsh.info/>`_. 
In the example the ``elementLabels`` are set to -1 on :math:`\Gamma^D` (blue in the below figure) and 1 on :math:`\Omega` (yellow). 
The plot is based on a domain saved in ``examples/docs/mesh.pkl``.

::

    import pickle
    import matplotlib.pyplot as plt

    mesh_arrays = pickle.load(open("mesh.pkl", "rb"))
    elements, vertices, elementLabels = mesh_arrays["elements"], mesh_arrays["vertices"], mesh_arrays["elementLabels"]

    plt.tripcolor(vertices[:, 0], vertices[:, 1], elements, facecolors=elementLabels, edgecolors="black", linewidth=.2, alpha=.7)
    plt.title(r"Domain $\Omega \cup \Gamma^D$")
    plt.gca().set_aspect('equal')
    plt.savefig("domain.png")
    plt.close()


.. image:: ../../examples/docs/domain.png
	:width: 400
	
Settings
------------

The ``kernel`` is a dictionary which contains the keys ``function``, ``horizon``, ``outputdim`` and possibly   ``fractional_s``. 
A list of implemented kernel functions can be found via ``show_options()``.
The ``horizon`` determines the interaction radius of the kernel and ``outputdim`` tells whether the kernel is scalar or matrix-valued. Given the kernel has a singularity which requires regularizing integral transformations we need to determine the parameter s of the singularity in ``fractional_s``. The quantity is required for the special integration routines ``fractional`` and ``weakFractional`` which need to be applied to those kernels.

**Example** We choose a scalar kernel of truncated fractional type

.. math:: 

  C(\mathbf{x}, \mathbf{y}) = 
  \frac{2-2s}{\pi\delta^{2-2s}}
  \frac{1}{|\mathbf{x}-\mathbf{y}|^{d + 2s}}
  

where :math:`d=2` and :math:`s=0.4`. The corresponding dictionary reads as

::

    kernel = {
        "function": "fractional",
        "horizon": 0.1,
        "outputdim": 1,
        "fractional_s": 0.4
    }

As optional parameter it is possible to hand over varying kernel coefficients with the key ``"Theta"``. Of course, the implemented kernel has to support the usage of the given information. 
This is so, for example, for the kernel functions ``"theta"`` and  ``"sparsetheta"`` which 
evaluate

.. math::
    \frac{4}{\pi \delta^{4}} \Theta(x,y)

and

.. math::
    \frac{4}{\pi \delta^{4}} (1 - \Theta(x,y)),

respectively. The coefficients ``Theta`` are expected to be of type ``scipy.sparse.csr_matrix``. However, there are no other restrictions, so that any use of the coefficients inside of the kernel function is possible.
Note also that new kernels can be defined.

The ``forcing``-term has a similar structure and contains a key ``function`` indicating the forcing term.
Again, a list of implemented kernel functions can be found via ``show_options()``. 
Note that the assembly of the right hand side is standard and nlfem offers this functionality for convenience only.

**Example** We choose :math:`f(x) = -2(x_1+1)` which can be selected by

::

     forcing = {"function": "linear"}

All other settings are stored in the dictionary ``conf``. 
The truncation routines are selected by the ``method`` given in ``approxBalls``. 
The quadrature rule which is used in the truncation routine is indicated by numpy arrays of quadrature points and weights. 
For kernels which exhibit a singularity it is possible to select another integration routine specifically for touching elements 
(e.g. ``fractional``, ``weakFractional``). However, the choice does not affect the quadrature for nonsingular kernels. 
The function ``show_options()`` prints a list of the singular kernels.

**Example** We discretize the problem with a continuous Galerkin ansatz space, and use a retriangulation with caps to approximate 
the truncation of the interaction set. We choose a seven-point quadrature rule with points ``Px`` and weights ``dx`` which is advantageous for the  approximation of the truncated interaction set of a finite element, as the points lie close to the edges of the element [DElia2021]_. 

::

	import numpy as np
	Px = np.array([[0.33333333333333, 0.33333333333333],
			[0.47014206410511, 0.47014206410511],
			[0.47014206410511, 0.05971587178977],
			[0.05971587178977, 0.47014206410511],
			[0.10128650732346, 0.10128650732346],
			[0.10128650732346, 0.79742698535309],
			[0.79742698535309, 0.10128650732346]])
	dx = 0.5 * np.array([0.22500000000000,
			0.13239415278851,
			0.13239415278851,
			0.13239415278851,
			0.12593918054483,
			0.12593918054483,
			0.12593918054483])

.. image:: ../../examples/docs/quadrature.png
	:width: 300

The example configuration is set by the dictionary:

::

   conf = {
        "ansatz": "CG",
        "approxBalls": {
            "method": "retriangulate"
        },
        "quadrature": {
            "outer": {
                "points": Px,
                "weights": dx
            },
            "inner": {
                "points": Px,
                "weights": dx
            },
            "touchingElements": {
            	 "method": "fractional",
            	 "ntensorGaussPoints": 4
            } 
        },
        "verbose": True
    }
    
**Note** All options for the settings in ``kernel``, ``function`` and ``conf`` can be printed by ``nlfem.show_options()``.
Empty dictionaries in the required form can be obtained from ``nlfem.get_empty_settings()``.

Assembly
------------
Given the above settings the assembly of the nonlocal stiffness matrix is effected via

::

    import nlfem
    mesh, A = nlfem.stiffnessMatrix_fromArray(elements, elementLabels, vertices, kernel, conf)
    f = nlfem.loadVector(mesh, forcing, conf)
	
The dictionary ``mesh`` contains information about the labeling of ``elements``, ``vertices`` and ``dof``. 
The labels of the degrees of freedom depend on the kernel (scalar or matrix-valued) and the ansatz space (continuous (CG) or discontinuous (DG) P1 elements). 
Therefore the function ``nlfem.stiffnessMatrix_fromArray()`` automatically deduces ``vertexLabels`` and ``dofLabels``. 
A vertex is labeled according to the labels of the elements it belongs to and it is given the smallest of those element labels.
The dof labels are identical to the vertexlabels for CG ansatz spaces and scalar-valued kernels. 
For DG ansatz spaces the labels are directly obtained from the element labels.

The labels of the degrees of freedom on the domain :math:`\Omega` and the Dirchlet domain :math:`\Gamma^D` are derived from the element labels and stored in ``mesh`` as ``dofLabels``. 

::

	Omega = mesh["dofLabels"] > 0
	Gamma = mesh["dofLabels"] < 0

**Example** In the given example we have a CG ansatz space and a scalar-valued kernel so that the ``dofLabels`` are identical to the ``vertexLabels``.

::

	plt.triplot(vertices[:, 0], vertices[:, 1], elements, color="black", linewidth=.2, alpha=.7)
	plt.scatter(vertices[:, 0], vertices[:, 1], c=mesh["dofLabels"])
	plt.title("DoF labels")
	plt.gca().set_aspect('equal')
	plt.savefig("doflabels.png")
	plt.close()

.. image:: ../../examples/docs/doflabels.png
	:width: 400
	
The solution of a Dirichlet-type problem requires the definition of Dirichlet data. The function below implements :math:`g(x) = x_1^2 x_2 + x_2^2`.

::

	def g(vertices):
		return vertices[:, 0] ** 2 * vertices[:, 1] + vertices[:, 1] ** 2


Ultimately we can solve the nonlocal Dirichlet-type problem 
via

::

    from scipy.sparse.linalg import cg
    f_O = f[Omega]
    A_OO = A[Omega][:, Omega]
    A_OD = A[Omega][:, Gamma]
    g_D = g(vertices[Gamma])
    f_O -= A_OD.dot(g_D)

    u = np.zeros(vertices.shape[0], dtype=float)
    u[Omega] = cg(A_OO, f_O)[0]
    u[Gamma] = g_D

::

	fig = plt.figure(figsize=plt.figaspect(0.5))
	ax = fig.add_subplot(1, 2, 1, projection='3d')
	surf = ax.plot_trisurf(vertices[:, 0],vertices[:, 1], u, cmap=plt.cm.Spectral)
	plt.savefig("solution.png")
	plt.close()
	
	
.. image:: ../../examples/docs/solution.png

Defining integral kernels
---------------------------

New kernels can be added to nlfem in the following way. Declare your kernel ``kernel_myKernel()`` in ``src/model.h`` and implement it in ``src/model.cpp``. You kernel function has to exactly meet the interface of the
function pointer ``model_kernel``.

:: 

	void (*model_kernel) (const double * x, long labelx,
	const double * y,long labely,
	const MeshType &mesh, double * kernel_val);

In order to be able to choose the kernel in the Python interface add an entry 

::

	{"myKernelName", kernel_myKernel}
	
to the map ``lookup_kernel`` in ``src/Cassemble.cpp``. 

.. [DElia2021] D’Elia, M., Gunzburger, M., & Vollmann, C. (2021). A cookbook for approximating Euclidean balls and for quadrature rules in finite element methods for nonlocal problems. Mathematical Models and Methods in Applied Sciences, 31(08), 1505-1567.

.. [Klar2022] Klar, M., Vollmann, C., & Schulz, V. (2022). nlfem: A flexible 2d Fem Code for Nonlocal Convection-Diffusion and Mechanics. arXiv preprint arXiv:2207.03921.

Python interface
==================

.. autofunction:: nlfem.show_options()
.. autofunction:: nlfem.get_empty_settings()
.. autofunction:: nlfem.stiffnessMatrix_fromArray
.. autofunction:: nlfem.stiffnessMatrix
.. autofunction:: nlfem.loadVector

Indices and tables
==================

* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`