Boundary conditions

class BoundaryCondition(object):
    """
    Base class for boundary conditions in a LBM simulation.

    This class provides a general structure for implementing boundary conditions. It includes methods for preparing the
    boundary attributes and for applying the boundary condition. Specific boundary conditions should be implemented as
    subclasses of this class, with the `apply` method overridden as necessary.

    Attributes
    ----------
    lattice (Lattice): The lattice used in the simulation.

    nx (int): The number of nodes in the x direction.

    ny (int): The number of nodes in the y direction.

    nz (int): The number of nodes in the z direction.

    dim (int): The number of dimensions in the simulation (2 or 3).

    precision_policy (PrecisionPolicy): The precision policy used in the simulation.

    indices (array-like): The indices of the boundary nodes.

    name (str or None): The name of the boundary condition. This should be set in subclasses.

    is_solid (bool): Whether the boundary condition is for a solid boundary. This should be set in subclasses.

    is_dynamic (bool): Whether the boundary condition is dynamic (changes over time). This should be set in subclasses.

    needs_extra_configuration (bool): Whether the boundary condition requires extra configuration. This should be set in subclasses.

    implementation_step (str): The lattice Boltzmann algorithm step at which the boundary condition is applied.
    """

    def __init__(self, indices, grid_info, precision_policy):
        self.lattice = grid_info["lattice"]
        self.nx = grid_info["nx"]
        self.ny = grid_info["ny"]
        self.nz = grid_info["nz"]
        self.dim = grid_info["dim"]
        self.precision_policy = precision_policy
        self.indices = indices
        self.name = None
        self.is_solid = False
        self.is_dynamic = False
        self.needs_extra_configuration = False
        self.implementation_step = "PostStreaming"

    def create_local_mask_and_normal_arrays(self, grid_mask):
        """
        Creates local mask and normal arrays for the boundary condition.

        Parameters
        ----------
        grid_mask (array-like): The grid mask for the lattice.

        Returns
        -------
        None

        Notes
        -----
        This method creates local mask and normal arrays for the boundary condition based on the grid mask.
        If the boundary condition requires extra configuration, the `configure` method is called.
        """

        if self.needs_extra_configuration:
            boundaryMask = self.get_boundary_mask(grid_mask)
            self.configure(boundaryMask)
            self.needs_extra_configuration = False

        boundaryMask = self.get_boundary_mask(grid_mask)
        self.normals = self.get_normals(boundaryMask)
        self.imissing, self.iknown = self.get_missing_indices(boundaryMask)
        self.imissing_mask, self.iknown_mask, self.imiddle_mask = self.get_missing_mask(boundaryMask)

        return

    def get_boundary_mask(self, grid_mask):
        """
        Add jax.device_count() to the self.indices in x-direction, and 1 to the self.indices other directions
        This is to make sure the boundary condition is applied to the correct nodes as grid_mask is
        expanded by (jax.device_count(), 1, 1)

        Parameters
        ----------
        grid_mask (array-like): The grid mask for the lattice.

        Returns
        -------
        boundaryMask : (array-like)
        """
        shifted_indices = np.array(self.indices)
        shifted_indices[0] += device_count()
        shifted_indices[1:] += 1
        # Convert back to tuple
        shifted_indices = tuple(shifted_indices)
        boundaryMask = np.array(grid_mask[shifted_indices])

        return boundaryMask

    def configure(self, boundary_mask):
        """
        Configures the boundary condition.

        Parameters
        ----------
        boundary_mask (array-like): The grid mask for the boundary voxels.

        Returns
        -------
        None

        Notes
        -----
        This method should be overridden in subclasses if the boundary condition requires extra configuration.
        """
        return

    @partial(jit, static_argnums=(0, 3), inline=True)
    def prepare_populations(self, fout, fin, implementation_step):
        """
        Prepares the distribution functions for the boundary condition.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The incoming distribution functions.

        fin (jax.numpy.ndarray): The outgoing distribution functions.

        implementation_step (str): The step in the lattice Boltzmann method algorithm at which the preparation is
            applied.

        Returns
        -------
        (jax.numpy.ndarray): The prepared distribution functions.

        Notes
        -----
        This method should be overridden in subclasses if the boundary condition requires preparation of the distribution functions during post-collision or post-streaming.
        See ExtrapolationBoundaryCondition for an example.
        """
        return fout

    def get_normals(self, boundary_mask):
        """
        Calculates the normal vectors at the boundary nodes.

        Parameters
        ----------
        boundary_mask (array-like): The boundary mask for the lattice.

        Returns
        -------
        (array-like): The normal vectors at the boundary nodes.

        Notes
        -----
        This method calculates the normal vectors by dotting the boundary mask with the main lattice directions.
        """
        main_c = self.lattice.c.T[self.lattice.main_indices]
        m = boundary_mask[..., self.lattice.main_indices]
        normals = -np.dot(m, main_c)
        return normals

    def get_missing_indices(self, boundary_mask):
        """
        Returns two int8 arrays the same shape as boundary_mask. The non-zero entries of these arrays indicate missing
        directions that require BCs (imissing) as well as their corresponding opposite directions (iknown).

        Parameters
        ----------
        boundary_mask (array-like): The boundary mask for the lattice.

        Returns
        -------
        (tuple of array-like): The missing and known indices for the boundary condition.

        Notes
        -----
        This method calculates the missing and known indices based on the boundary mask. The missing indices are the
        non-zero entries of the boundary mask, and the known indices are their corresponding opposite directions.
        """

        # Find imissing, iknown 1-to-1 corresponding indices
        # Note: the "zero" index is used as default value here and won't affect BC computations
        nbd = len(self.indices[0])
        imissing = np.vstack([np.arange(self.lattice.q, dtype="uint8")] * nbd)
        iknown = np.vstack([self.lattice.opp_indices] * nbd)
        imissing[~boundary_mask] = 0
        iknown[~boundary_mask] = 0
        return imissing, iknown

    def get_missing_mask(self, boundary_mask):
        """
        Returns three boolean arrays the same shape as boundary_mask.
        Note: these boundary masks are useful for reduction (eg. summation) operators of selected q-directions.

        Parameters
        ----------
        boundary_mask (array-like): The boundary mask for the lattice.

        Returns
        -------
        (tuple of array-like): The missing, known, and middle masks for the boundary condition.

        Notes
        -----
        This method calculates the missing, known, and middle masks based on the boundary mask. The missing mask
        is the boundary mask, the known mask is the opposite directions of the missing mask, and the middle mask
        is the directions that are neither missing nor known.
        """
        # Find masks for imissing, iknown and imiddle
        imissingMask = boundary_mask
        iknownMask = imissingMask[:, self.lattice.opp_indices]
        imiddleMask = ~(imissingMask | iknownMask)
        return imissingMask, iknownMask, imiddleMask

    @partial(jit, static_argnums=(0,))
    def apply(self, fout, fin):
        """
        Applies the boundary condition.

        Parameters
        ----------
        fout (jax.Array): Output distribution functions.

        fin (jax.Array): Input distribution functions.

        Returns
        -------
        None

        Notes
        -----
        This method should be overridden in subclasses to implement the specific boundary condition. The method should
        modify the output distribution functions in place to apply the boundary condition.
        """
        pass

    @partial(jit, static_argnums=(0,))
    def equilibrium(self, rho, u):
        """
        Compute equilibrium distribution function.

        Parameters
        ----------
        rho (jax.numpy.ndarray): The density at each node in the lattice.

        u (jax.numpy.ndarray): The velocity at each node in the lattice.

        Returns
        -------
        (jax.numpy.ndarray): The equilibrium distribution function at each node in the lattice.

        Notes
        -----
        This method computes the equilibrium distribution function based on the density and velocity. The computation is
        performed in the compute precision specified by the precision policy. The result is not cast to the output precision as
        this is function is used inside other functions that require the compute precision.
        """
        rho, u = self.precision_policy.cast_to_compute((rho, u))
        c = jnp.array(self.lattice.c, dtype=self.precision_policy.compute_dtype)
        cu = 3.0 * jnp.dot(u, c)
        usqr = 1.5 * jnp.sum(u**2, axis=-1, keepdims=True)
        feq = rho * self.lattice.w * (1.0 + 1.0 * cu + 0.5 * cu**2 - usqr)

        return feq

    @partial(jit, static_argnums=(0,))
    def momentum_flux(self, fneq):
        """
        Compute the momentum flux.

        Parameters
        ----------
        fneq (jax.numpy.ndarray): The non-equilibrium distribution function at each node in the lattice.

        Returns
        -------
        (jax.numpy.ndarray): The momentum flux at each node in the lattice.

        Notes
        -----
        This method computes the momentum flux by dotting the non-equilibrium distribution function with the lattice
        direction vectors.
        """
        return jnp.dot(fneq, self.lattice.cc)

    @partial(jit, static_argnums=(0,))
    def momentum_exchange_force(self, f_poststreaming, f_postcollision):
        """
        Using the momentum exchange method to compute the boundary force vector exerted on the solid geometry
        based on [1] as described in [3]. Ref [2] shows how [1] is applicable to curved geometries only by using a
        bounce-back method (e.g. Bouzidi) that accounts for curved boundaries.
        NOTE: this function should be called after BC's are imposed.
        [1] A.J.C. Ladd, Numerical simulations of particular suspensions via a discretized Boltzmann equation.
            Part 2 (numerical results), J. Fluid Mech. 271 (1994) 311-339.
        [2] R. Mei, D. Yu, W. Shyy, L.-S. Luo, Force evaluation in the lattice Boltzmann method involving
            curved geometry, Phys. Rev. E 65 (2002) 041203.
        [3] Caiazzo, A., & Junk, M. (2008). Boundary forces in lattice Boltzmann: Analysis of momentum exchange
            algorithm. Computers & Mathematics with Applications, 55(7), 1415-1423.

        Parameters
        ----------
        f_poststreaming (jax.numpy.ndarray): The post-streaming distribution function at each node in the lattice.

        f_postcollision (jax.numpy.ndarray): The post-collision distribution function at each node in the lattice.

        Returns
        -------
        (jax.numpy.ndarray): The force exerted on the solid geometry at each boundary node.

        Notes
        -----
        This method computes the force exerted on the solid geometry at each boundary node using the momentum exchange method.
        The force is computed based on the post-streaming and post-collision distribution functions. This method
        should be called after the boundary conditions are imposed.
        """
        c = jnp.array(self.lattice.c, dtype=self.precision_policy.compute_dtype)
        nbd = len(self.indices[0])
        bindex = np.arange(nbd)[:, None]
        phi = f_postcollision[self.indices][bindex, self.iknown] + f_poststreaming[self.indices][bindex, self.imissing]
        force = jnp.sum(c[:, self.iknown] * phi, axis=-1).T
        return force

class BounceBack(BoundaryCondition):
    """
    Bounce-back boundary condition for a lattice Boltzmann method simulation.

    This class implements a full-way bounce-back boundary condition, where particles hitting the boundary are reflected
    back in the direction they came from. The boundary condition is applied after the collision step.

    Attributes
    ----------
    name (str): The name of the boundary condition. For this class, it is "BounceBackFullway".

    implementation_step (str): The step in the lattice Boltzmann method algorithm at which the boundary condition is applied. For this class, it is "PostCollision".

    theta (jax.numpy.ndarray; Default: None): Contact angle, applied for multiphase flows and only set for wall boundary conditions.

    phi (jax.numpy.ndarray; Default: None): Contact angle parameter phi, applied for multiphase flows and only set for wall boundary conditions.

    delta_rho (jax.numpy.ndarray; Default: None): Contact angle parameter delta_rho, applied for multiphase flows and only set for wall boundary conditions.
    """

    def __init__(self, indices, grid_info, precision_policy, theta=None, phi=None, delta_rho=None):
        super().__init__(indices, grid_info, precision_policy)
        self.name = "BounceBackFullway"
        self.implementation_step = "PostCollision"
        self.theta = theta
        self.phi = phi
        self.delta_rho = delta_rho

    @partial(jit, static_argnums=(0,))
    def apply(self, fout, fin):
        """
        Applies the bounce-back boundary condition.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The output distribution functions.

        fin (jax.numpy.ndarray): The input distribution functions.

        Returns
        -------
        (jax.numpy.ndarray): The modified output distribution functions after applying the boundary condition.

        Notes
        -----
        This method applies the bounce-back boundary condition by reflecting the input distribution functions at the
        boundary nodes in the opposite direction.
        """

        return fin[self.indices][..., self.lattice.opp_indices]

class BounceBackMoving(BoundaryCondition):
    """
    Moving bounce-back boundary condition for a lattice Boltzmann method simulation.

    This class implements a moving bounce-back boundary condition, where particles hitting the boundary are reflected
    back in the direction they came from, with an additional velocity due to the movement of the boundary. The boundary
    condition is applied after the collision step.

    Attributes
    ----------
    name (str): The name of the boundary condition. For this class, it is "BounceBackFullwayMoving".

    implementation_step (str): The step in the lattice Boltzmann method algorithm at which the boundary condition is applied. For this class, it is "PostCollision".

    is_dynamic (bool): Whether the boundary condition is dynamic (changes over time). For this class, it is True.

    update_function (function): A function that updates the boundary condition. For this class, it is a function that updates the boundary
    condition based on the current time step. The signature of the function is `update_function(time) -> (indices, vel)`,

    theta (jax.numpy.ndarray; Default: None): Contact angle, applied for multiphase flows and only set for wall boundary conditions.

    phi (jax.numpy.ndarray; Default: None): Contact angle parameter phi, applied for multiphase flows and only set for wall boundary conditions.

    delta_rho (pytree of jax.numpy.ndarray; Default: None): Contact angle parameter delta_rho, applied for multiphase flows and only set for wall boundary conditions.
    """

    def __init__(self, grid_info, precision_policy, update_function=None, theta=None, phi=None, delta_rho=None):
        # We get the indices at time zero to pass to the parent class for initialization
        indices, _ = update_function(0)
        super().__init__(indices, grid_info, precision_policy)
        self.name = "BounceBackFullwayMoving"
        self.implementation_step = "PostCollision"
        self.is_dynamic = True
        self.update_function = jit(update_function)
        self.theta = theta
        self.phi = phi
        self.delta_rho = delta_rho

    @partial(jit, static_argnums=(0,))
    def apply(self, fout, fin, time):
        """
        Applies the moving bounce-back boundary condition.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The output distribution functions.

        fin (jax.numpy.ndarray): The input distribution functions.

        time (int): The current time step.

        Returns
        -------
        (jax.numpy.ndarray): The modified output distribution functions after applying the boundary condition.
        """
        indices, vel = self.update_function(time)
        c = jnp.array(self.lattice.c, dtype=self.precision_policy.compute_dtype)
        cu = 6.0 * self.lattice.w * jnp.dot(vel, c)
        return fout.at[indices].set(fin[indices][..., self.lattice.opp_indices] - cu)

class BounceBackHalfway(BoundaryCondition):
    """
    Halfway bounce-back boundary condition for a lattice Boltzmann method simulation.

    This class implements a halfway bounce-back boundary condition. The boundary condition is applied after
    the streaming step.

    Attributes
    ----------
    name (str): The name of the boundary condition. For this class, it is "BounceBackHalfway".

    implementation_step (str): The step in the lattice Boltzmann method algorithm at which the boundary condition is applied. For this class, it is "PostStreaming".

    needs_extra_configuration (bool): Whether the boundary condition needs extra configuration before it can be applied. For this class, it is True.

    is_solid (bool): Whether the boundary condition represents a solid boundary. For this class, it is True.

    solid_indices (tuple): Original solid-node indices, stored by configure before self.indices is shifted to the
    adjacent fluid nodes. Used by multiphase wetting schemes.

    vel (array-like): The prescribed value of velocity vector for the boundary condition. No-slip BC is assumed if vel=None (default).

    theta (jax.numpy.ndarray; Default: None): Contact angle, applied for multiphase flows and only set for wall boundary conditions.

    phi (jax.numpy.ndarray; Default: None): Contact angle parameter phi, applied for multiphase flows and only set for wall boundary conditions.

    delta_rho (pytree of jax.numpy.ndarray; Default: None): Contact angle parameter delta_rho, applied for multiphase flows and only set for wall boundary conditions.
    """

    def __init__(self, indices, grid_info, precision_policy, vel=None, theta=None, phi=None, delta_rho=None):
        super().__init__(indices, grid_info, precision_policy)
        self.name = "BounceBackHalfway"
        self.implementation_step = "PostStreaming"
        self.needs_extra_configuration = True
        self.is_solid = True
        self.vel = vel
        self.theta = theta
        self.phi = phi
        self.delta_rho = delta_rho

    def configure(self, boundary_mask):
        """
        Configures the boundary condition.

        Parameters
        ----------
        boundary_mask (array-like): The grid mask for the boundary voxels.

        Returns
        -------
        None

        Notes
        -----
        This method performs an index shift for the halfway bounce-back boundary condition. It updates the indices of
        the boundary nodes to be the indices of fluid nodes adjacent of the solid nodes.
        """
        # Keep the original solid-node indices; multiphase wetting schemes need them
        # after self.indices is shifted to the adjacent fluid nodes below.
        self.solid_indices = self.indices
        # Perform index shift for halfway BB.
        hasFluidNeighbour = ~boundary_mask[:, self.lattice.opp_indices]
        nbd_orig = len(self.indices[0])
        idx = np.array(self.indices).T
        idx_trg = []
        for i in range(self.lattice.q):
            idx_trg.append(idx[hasFluidNeighbour[:, i], :] + self.lattice.c[:, i])
        indices_new = np.unique(np.vstack(idx_trg), axis=0)
        self.indices = tuple(indices_new.T)
        nbd_modified = len(self.indices[0])
        if (nbd_orig != nbd_modified) and self.vel is not None:
            vel_avg = np.mean(self.vel, axis=0)
            self.vel = jnp.zeros(indices_new.shape, dtype=self.precision_policy.compute_dtype) + vel_avg
            warnings.warn(
                "A constant averaged velocity vector is imposed at all boundary cells.",
                UserWarning,
                stacklevel=2,
            )

        return

    @partial(jit, static_argnums=(0, 3), inline=True)
    def prepare_populations(self, fout, fin, implementation_step):
        """
        Pin solid-node populations to the rest equilibrium (rho = 1, u = 0).

        Halfway bounce-back never constrains the populations at the solid nodes, so they evolve freely and can diverge.
        This is harmless in single phase simulations but multiphase wetting schemes read wall densities, so the solid
        nodes are reset after every streaming step. The wetted wall density used in the force computation is overwritten
        by apply_contact_angle, so the pinned value never enters the wetting force directly.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The post-streaming or post-collision distribution functions.

        fin (jax.numpy.ndarray): The pre-collision or post-collision distribution functions.

        implementation_step (str): The step at which the preparation is applied.

        Returns
        -------
        (jax.numpy.ndarray): The distribution functions with solid nodes reset after streaming.
        """
        if implementation_step == "PostStreaming":
            return fout.at[self.solid_indices].set(self.precision_policy.cast_to_output(self.lattice.w))
        return fout

    @staticmethod
    def reflect_missing(fbd, fin_bd, imissing, iknown):
        """
        Core halfway bounce-back math: set each missing direction to the post-streaming value coming from its
        known (opposite) direction. Pure function of its arguments (no BC-instance state), so it is reusable
        both by apply() (global indices) and by a local (per-shard) bounce-back kernel operating on a gathered
        block with the same (n, q) shape.

        Parameters
        ----------
        fbd (jax.numpy.ndarray): Post-streaming distribution functions at the boundary nodes, shape (n, q).

        fin_bd (jax.numpy.ndarray): Post-streaming distribution functions at the same boundary nodes, shape (n, q).

        imissing (array-like): Missing-direction indices, shape (n, q).

        iknown (array-like): Known (opposite) direction indices, shape (n, q).

        Returns
        -------
        (jax.numpy.ndarray): fbd with each missing direction set from its known direction.
        """
        bindex = jnp.arange(fbd.shape[0])[:, None]
        return fbd.at[bindex, imissing].set(fin_bd[bindex, iknown])

    @staticmethod
    def velocity_correction(fbd, imissing, iknown, vel, w, c):
        """
        Core halfway bounce-back velocity-forcing math, factored out of impose_boundary_vel so it is reusable by
        a local (per-shard) bounce-back kernel. Pure function of its arguments; vel=0 is a safe no-op, so a
        local kernel merging several boundary conditions can always call this instead of conditionally skipping
        it per boundary condition.

        Parameters
        ----------
        fbd (jax.numpy.ndarray): Distribution functions at the boundary nodes, shape (n, q).

        imissing (array-like): Missing-direction indices, shape (n, q).

        iknown (array-like): Known (opposite) direction indices, shape (n, q).

        vel (jax.numpy.ndarray): Prescribed velocity vector, shape (n, dim).

        w (jax.numpy.ndarray): Lattice weights, shape (q,).

        c (jax.numpy.ndarray): Lattice velocity vectors, shape (dim, q).

        Returns
        -------
        (jax.numpy.ndarray): fbd with the velocity correction applied.
        """
        bindex = jnp.arange(fbd.shape[0])[:, None]
        cu = 6.0 * w * jnp.dot(vel, c)
        return fbd.at[bindex, imissing].add(-cu[bindex, iknown])

    @partial(jit, static_argnums=(0,))
    def impose_boundary_vel(self, fbd, bindex):
        c = jnp.array(self.lattice.c, dtype=self.precision_policy.compute_dtype)
        return self.velocity_correction(fbd, self.imissing, self.iknown, self.vel, self.lattice.w, c)

    @partial(jit, static_argnums=(0,))
    def apply(self, fout, fin):
        """
        Applies the halfway bounce-back boundary condition.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The output distribution functions.

        fin (jax.numpy.ndarray): The input distribution functions.

        Returns
        -------
        (jax.numpy.ndarray): The modified output distribution functions after applying the boundary condition.
        """
        fbd = self.reflect_missing(fout[self.indices], fin[self.indices], self.imissing, self.iknown)
        if self.vel is not None:
            bindex = np.arange(len(self.indices[0]))[:, None]
            fbd = self.impose_boundary_vel(fbd, bindex)
        return fbd

class EquilibriumBC(BoundaryCondition):
    """
    Equilibrium boundary condition for a lattice Boltzmann method simulation.

    This class implements an equilibrium boundary condition, where the distribution function at the boundary nodes is
    set to the equilibrium distribution function. The boundary condition is applied after the streaming step.

    Attributes
    ----------
    name (str): The name of the boundary condition. For this class, it is "EquilibriumBC".

    implementation_step (str): The step in the lattice Boltzmann method algorithm at which the boundary condition is applied. For this class, it is "PostStreaming".

    out (jax.numpy.ndarray): The equilibrium distribution function at the boundary nodes.
    """

    def __init__(self, indices, grid_info, precision_policy, rho, u):
        super().__init__(indices, grid_info, precision_policy)
        self.out = self.precision_policy.cast_to_output(self.equilibrium(rho, u))
        self.name = "EquilibriumBC"
        self.implementation_step = "PostStreaming"

    @partial(jit, static_argnums=(0,))
    def apply(self, fout, fin):
        """
        Applies the equilibrium boundary condition.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The output distribution functions.

        fin (jax.numpy.ndarray): The input distribution functions.

        Returns
        -------
        (jax.numpy.ndarray): The modified output distribution functions after applying the boundary condition.

        Notes
        -----
        This method applies the equilibrium boundary condition by setting the output distribution functions at the
        boundary nodes to the equilibrium distribution function.
        """
        return self.out

class DoNothing(BoundaryCondition):
    """
    Preserve post-collision populations at selected boundary nodes.

    Parameters
    ----------
    indices (tuple of numpy.ndarray): Boundary-node indices.

    grid_info (dict): Grid and lattice metadata.

    precision_policy (PrecisionPolicy): Compute and output precision policy.

    Notes
    -----
    Streaming is skipped at these nodes by returning the post-collision
    populations after streaming. This avoids values wrapping into the domain
    from the opposite side of the rolled population array.
    """

    def __init__(self, indices, grid_info, precision_policy):
        super().__init__(indices, grid_info, precision_policy)
        self.name = "DoNothing"
        self.implementation_step = "PostStreaming"

    @partial(jit, static_argnums=(0,))
    def apply(self, fout, fin):
        """
        Applies the do-nothing boundary condition.

        Parameters
        ----------
        fout (jax.numpy.ndarray): The output distribution functions.

        fin (jax.numpy.ndarray): The input distribution functions.

        Returns
        -------
        jax.numpy.ndarray
            The modified output distribution functions after applying the boundary condition.

        Notes
        -----
        This method applies the do-nothing boundary condition by simply returning the input distribution functions at the
        boundary nodes.
        """
        return fin[self.indices]