.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/example_nonlocal.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_example_nonlocal.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_example_nonlocal.py:


Nonlocal problems
=================

.. GENERATED FROM PYTHON SOURCE LINES 15-50

I this second example, we will assemble and solve several nonlocal equations.
The full code of this example can be found in `examples/example_nonlocal.py <https://github.com/sandialabs/PyNucleus/blob/master/examples/example_nonlocal.py>`_ in the PyNucleus repository.

PyNucleus can assemble operators of the form

.. math::

   \mathcal{L}u(x) = \operatorname{p.v.} \int_{\mathbb{R}^d} [u(y)-u(x)] \gamma(x, y) dy.

The kernel :math:`\gamma` is of the form

.. math::

   \gamma(x,y) = \frac{\phi(x, y)}{|x-y|^{\beta(x,y)}} \chi_{\mathcal{N}(x)}(y).

Here, :math:`\phi` is a positive function, and :math:`\chi` is the indicator function.
The interaction neighborhood :math:`\mathcal{N}(x) \subset \mathbb{R}^d` is often given as parametrized by the so-called horizon :math:`0<\delta\le\infty`:

.. math::

   \mathcal{N}(x) = \{y \in \mathbb{R}^d \text{ such that } ||x-y||_p<\delta\}.

where :math:`p \in \{1,2,\infty\}`. Other types of neighborhoods are also possible.

The singularity :math:`\beta` of the kernel crucially influences properties of the kernel:

- fractional type: :math:`\beta(x,y)=d+2s(x,y)`, where :math:`d` is the spatial dimension and :math:`s(x,y)` is the fractional order.
- constant type :math:`\beta(x,y)=0`
- peridynamic type :math:`\beta(x,y)=1`


A fractional kernel
-------------------

We start off by creating a fractional kernel with infinite horizon and constant fractional order :math:`s=0.75`.

.. GENERATED FROM PYTHON SOURCE LINES 50-56

.. code-block:: Python


    import matplotlib.pyplot as plt
    from time import time
    from PyNucleus import kernelFactory
    print(kernelFactory)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Available:
    'fractional', signature: '(dim, s, horizon=None, interaction=None, scaling=None, normalized=True, piecewise=True, phi=None, boundary=False, derivative=0, tempered=0.0, max_horizon=nan, manifold=False)'
    'indicator' with aliases ['constant'], signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'
    'inverseDistance' with aliases ['peridynamic', 'inverseOfDistance'], signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'
    'gaussian', signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'
    'exponential', signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'
    'polynomial', signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'
    'logInverseDistance', signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'
    'monomial', signature: '(dim, kernel, horizon, scaling=None, interaction=None, normalized=True, piecewise=True, phi=None, boundary=False, monomialPower=nan, variance=1.0, exponentialRate=1.0, a=1.0, max_horizon=nan)'





.. GENERATED FROM PYTHON SOURCE LINES 57-63

.. code-block:: Python


    from numpy import inf

    kernelFracInf = kernelFactory('fractional', dim=2, s=0.75, horizon=inf)
    print(kernelFracInf)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    kernel(fractional, 0.75, R^2, constantFractionalLaplacianScaling(0.75,inf -> 0.08558356484527617))




.. GENERATED FROM PYTHON SOURCE LINES 64-68

.. code-block:: Python


    plt.figure().gca().set_title('Fractional kernel')
    kernelFracInf.plot()




.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_001.png
   :alt: Fractional kernel
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 69-88

By default, kernels are normalized so that their local limits recover the classical Laplacian. This can be disabled by passing ``normalized=False``.


Nonlocal assembly
-----------------

We will be solving the problem

 .. math::

    -\mathcal{L} u &= f && \text{ in } \Omega=B(0,1)\subset\mathbb{R}^2, \\
    u &= 0 && \text{ in } \mathbb{R}^2 \setminus \Omega,

for constant forcing function :math:`f=1`.

First, we generate a mesh.
Instead of the ``meshFactory`` used in the previous example, we now use the ``nonlocalMeshFactory``.
The advantage is that this factory can generate meshes with appropriate interaction domains.
For this particular example, the factory will not generate any interaction domain, since the homogeneous Dirichlet condition on :math:`\mathbb{R}^2\setminus\Omega` can be enforced via a boundary integral.

.. GENERATED FROM PYTHON SOURCE LINES 88-96

.. code-block:: Python


    from PyNucleus import nonlocalMeshFactory, HOMOGENEOUS_DIRICHLET

    # Get a mesh that is appropriate for the problem, i.e. with the required interaction domain.
    meshFracInf, _ = nonlocalMeshFactory('disc', kernel=kernelFracInf, boundaryCondition=HOMOGENEOUS_DIRICHLET, hTarget=0.15)

    print(meshFracInf)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    mesh2d with 817 vertices and 1,536 cells




.. GENERATED FROM PYTHON SOURCE LINES 97-100

.. code-block:: Python

    plt.figure().gca().set_title('Mesh for fractional kernel')
    meshFracInf.plot()




.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_002.png
   :alt: Mesh for fractional kernel
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 101-105

Next, we obtain a piecewise linear, continuous DoFMap on the mesh, assemble the RHS and interpolate the known analytic solution.
We assemble the nonlocal operator by passing the kernel to the `assembleNonlocal` method of the DoFMap object.
The optional parameter `matrixFormat` determines what kind of linear operator is assembled.
We time the assembly of the operator as a dense matrix, and as a hierarchical matrix, and inspect the resulting objects.

.. GENERATED FROM PYTHON SOURCE LINES 105-117

.. code-block:: Python


    from PyNucleus import dofmapFactory, functionFactory

    dmFracInf = dofmapFactory('P1', meshFracInf)

    rhs = functionFactory('constant', 1.)
    exact_solution = functionFactory('solFractional', dim=2, s=0.75)

    b = dmFracInf.assembleRHS(rhs)
    u_exact = dmFracInf.interpolate(exact_solution)
    u = dmFracInf.zeros()








.. GENERATED FROM PYTHON SOURCE LINES 118-119

We assemble the operator in dense format.

.. GENERATED FROM PYTHON SOURCE LINES 119-126

.. code-block:: Python

    start = time()
    A_fracInf = dmFracInf.assembleNonlocal(kernelFracInf, matrixFormat='dense')

    print('Dense assembly took {}s'.format(time()-start))
    print(A_fracInf)
    print("Memory size: {} KB".format(A_fracInf.getMemorySize() >> 10))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Dense assembly took 2.5371222496032715s
    <721x721 Dense_LinearOperator>
    Memory size: 4061 KB




.. GENERATED FROM PYTHON SOURCE LINES 127-128

Then we assemble the operator in hierarchical matrix format.

.. GENERATED FROM PYTHON SOURCE LINES 128-136

.. code-block:: Python


    start = time()
    A_fracInf_h2 = dmFracInf.assembleNonlocal(kernelFracInf, matrixFormat='h2')

    print('Hierarchical assembly took {}s'.format(time()-start))
    print(A_fracInf_h2)
    print("Memory size: {} KB".format(A_fracInf_h2.getMemorySize() >> 10))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Hierarchical assembly took 2.611849546432495s
    <721x721 H2Matrix 0.362394 fill from near field, 0.104057 fill from tree, 0.167501 fill from clusters, 59 tree nodes, 164 far-field cluster pairs>
    Memory size: 2574 KB




.. GENERATED FROM PYTHON SOURCE LINES 137-143

It can be observed that both assembly routines take roughly the same amount of time.
The reason for this is that the operator itself has quite small dimensions.
For larger number of unknowns, we expect the hierarchical assembly scale like :math:`\mathcal{O}(N \log^{2d} N)`, whereas the dense assembly will scale at best like :math:`\mathcal{O}(N^2)`.
The memory usage of the hierarchical matrix is slightly better and scales similar to the complexity of assembly.

Similar to the local PDE example, we can then solve the resulting linear equation and compute the error in energy norm.

.. GENERATED FROM PYTHON SOURCE LINES 143-154

.. code-block:: Python


    from PyNucleus import solverFactory
    from numpy import sqrt

    solver = solverFactory('lu', A=A_fracInf, setup=True)
    solver(b, u)

    Hs_err = sqrt(abs(b.inner(u-u_exact)))

    print('Hs error: {}'.format(Hs_err))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Hs error: 0.058095118511890205




.. GENERATED FROM PYTHON SOURCE LINES 155-160

.. code-block:: Python


    plt.figure().gca().set_title('Numerical solution, fractional kernel')
    u.plot()





.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_003.png
   :alt: example nonlocal
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes3D: >



.. GENERATED FROM PYTHON SOURCE LINES 161-178

A finite horizon case with Dirichlet volume condition
-----------------------------------------------------

Next, we solve a nonlocal Poisson problem involving a constant kernel with finite horizon.
We will choose :math:`\gamma(x,y) \sim \chi_{\mathcal{N}(x)}(y)` for a neighborhood defined by the 2-norm and horizon :math:`\delta=0.2`, and solve

.. math::

   -\mathcal{L} u &= f && \text{ in } \Omega=[0,1]^2, \\
   u &= -(x_1^2 + x_2^2)/4 && \text{ in } \mathcal{I},

where the interaction domain is given by

.. math::

   \mathcal{I} := \{y\in\mathbb{R}^2\setminus\Omega \text{ such that } \exists x\in\Omega: \gamma(x,y)\neq 0\}.


.. GENERATED FROM PYTHON SOURCE LINES 178-182

.. code-block:: Python


    kernelConst = kernelFactory('constant', dim=2, horizon=0.2)
    print(kernelConst)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Kernel(indicator, |x-y|_2 <= 0.2, constantIntegrableScaling(0.2 -> 795.7747154594766))




.. GENERATED FROM PYTHON SOURCE LINES 183-187

.. code-block:: Python


    plt.figure().gca().set_title('Constant kernel')
    kernelConst.plot()




.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_004.png
   :alt: Constant kernel
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 188-196

.. code-block:: Python


    from PyNucleus import DIRICHLET

    meshConst, nIConst = nonlocalMeshFactory('square', kernel=kernelConst, boundaryCondition=DIRICHLET, hTarget=0.18)

    print(meshConst)
    print(nIConst)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    mesh2d with 569 vertices and 1,056 cells
    {'domain': squareIndicator, 'boundary': 1.0*1.0*1.0+-1.0*squareIndicator+-1.0*1.0*1.0+-1.0*squareIndicator, 'interaction': 1.0*1.0+-1.0*squareIndicator, 'tag': -128, 'zeroExterior': False}




.. GENERATED FROM PYTHON SOURCE LINES 197-198

The dictionary ``nIConst`` contains several indicator functions that we will use to distinguish between interior and interaction domain.

.. GENERATED FROM PYTHON SOURCE LINES 198-202

.. code-block:: Python


    plt.figure().gca().set_title('Mesh for constant kernel')
    meshConst.plot()




.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_005.png
   :alt: Mesh for constant kernel
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 203-206

We can observe that the ``nonlocalMeshFactory`` generated a mesh that includes the interaction domain.

We set up 2 DoFMaps this time, one for the unknown interior degrees of freedom and one for the Dirichlet volume condition.

.. GENERATED FROM PYTHON SOURCE LINES 206-213

.. code-block:: Python


    dmConst = dofmapFactory('P1', meshConst, nIConst['domain'])
    dmConstInteraction = dmConst.getComplementDoFMap()

    print(dmConst)
    print(dmConstInteraction)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    P1 DoFMap with 225 DoFs and 344 boundary DoFs.
    P1 DoFMap with 344 DoFs and 225 boundary DoFs.




.. GENERATED FROM PYTHON SOURCE LINES 214-216

We can see that they are complementary to each other.
Next, we assemble two matrices.

.. GENERATED FROM PYTHON SOURCE LINES 216-237

.. code-block:: Python


    A_const = dmConst.assembleNonlocal(kernelConst, matrixFormat='sparse')
    B_const = dmConst.assembleNonlocal(kernelConst, dm2=dmConstInteraction, matrixFormat='sparse')

    g = functionFactory('Lambda', lambda x: -(x[0]**2 + x[1]**2)/4)
    g_interp = dmConstInteraction.interpolate(g)

    b = dmConst.assembleRHS(rhs)-(B_const*g_interp)
    u = dmConst.zeros()

    solver = solverFactory('cg', A=A_const, setup=True)
    solver.maxIter = 1000
    solver.tolerance = 1e-8

    solver(b, u)

    u_global = u.augmentWithBoundaryData(g_interp)

    plt.figure().gca().set_title('Numerical solution, constant kernel')
    u_global.plot()




.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_006.png
   :alt: example nonlocal
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_006.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes3D: >



.. GENERATED FROM PYTHON SOURCE LINES 238-242

.. code-block:: Python


    plt.figure().gca().set_title('Analytic solution, constant kernel')
    u_global.dm.interpolate(g).plot()




.. image-sg:: /auto_examples/images/sphx_glr_example_nonlocal_007.png
   :alt: example nonlocal
   :srcset: /auto_examples/images/sphx_glr_example_nonlocal_007.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes3D: >



.. GENERATED FROM PYTHON SOURCE LINES 243-246

.. code-block:: Python


    print(A_const)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    <225x225 SSS_LinearOperator with 4437 stored elements>




.. GENERATED FROM PYTHON SOURCE LINES 247-250

.. code-block:: Python


    plt.show()
    # sphinx_gallery_thumbnail_number = 3








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 11.468 seconds)


.. _sphx_glr_download_auto_examples_example_nonlocal.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_nonlocal.ipynb <example_nonlocal.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_nonlocal.py <example_nonlocal.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_nonlocal.zip <example_nonlocal.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_