========================================
Demo: Using the Functions Module 
========================================

This demo illustrates the usage of the `functions` module in our library. Specifically, we will demonstrate 
how to create indicator functions for fictitious domains. 

First, we create two domains using the Mesh object. The first mesh, `domain` represents the entire computational 
domain, while the second mesh, `fictitious`, represents a subdomain within the main domain. 

.. math:: 
    
    \Omega_{fictitious} \subset \Omega_{domain}

.. code-block:: python

    import dolfinx

    from flatiron_tk.functions import build_field_scalar_function
    from flatiron_tk.functions import build_rank_indicator_function
    from flatiron_tk.mesh import RectMesh
    from mpi4py import MPI

    # Define the domain and fictitious region
    domain = RectMesh(0.0, 0.0, 10.0, 10.0, 1/20)
    fictitious = RectMesh(2.0, 2.0, 8.0, 8.0, 1/20)

Field Scalar (Indicator) Function
----------------------------------------

We will then create a function that indicates the presence of the fictitious domain within the main domain. 
The function will return 1 inside the fictitious domain and 0 outside. Note, the interior/exterior values can be
customized.

.. code-block:: python

    # Build a scalar function that is 1 inside the fictitious region and 0 outside
    inside_value = 1.0
    outside_value = 0.0
    scalar_function = build_field_scalar_function(domain, fictitious, inside_value, outside_value)

We then save the function to a file for visualization.

.. code-block:: python 

    # Write the scalar function to an XDMF file
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, 'scalar_function.xdmf', 'w') as xdmf:
        xdmf.write_mesh(domain.msh)
        xdmf.write_function(scalar_function)


Rank Indicator Function
----------------------------------------
Next, we create a rank indicator function that assigns a unique value to each MPI process. This is useful for
visualizing how the domain is partitioned across different processes while debugging parallel computations.

.. code-block:: python

    # Build a rank indicator function
    rank_function = build_rank_indicator_function(domain)

    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, 'rank_indicator_function.xdmf', 'w') as xdmf:
        xdmf.write_mesh(domain.msh)
        xdmf.write_function(rank_indicator_function)

Full Script
----------------------------------------
.. code-block:: python

    import dolfinx

    from flatiron_tk.functions import build_field_scalar_function
    from flatiron_tk.functions import build_rank_indicator_function
    from flatiron_tk.mesh import RectMesh
    from mpi4py import MPI

    # Define the domain and fictitious region
    domain = RectMesh(0.0, 0.0, 10.0, 10.0, 1/20)
    fictitious = RectMesh(2.0, 2.0, 8.0, 8.0, 1/20)

    # Build a scalar function that is 1 inside the fictitious region and 0 outside
    inside_value = 1.0
    outside_value = 0.0
    scalar_function = build_field_scalar_function(domain, fictitious, inside_value, outside_value)

    # Write the scalar function to an XDMF file
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, 'scalar_function.xdmf', 'w') as xdmf:
        xdmf.write_mesh(domain.msh)
        xdmf.write_function(scalar_function)

    # This function builds a rank indicator function for the mesh 
    rank_indicator_function = build_rank_indicator_function(domain)
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, 'rank_indicator_function.xdmf', 'w') as xdmf:
        xdmf.write_mesh(domain.msh)
        xdmf.write_function(rank_indicator_function)