Cumulative integral of an interpolated expression

dolfinx
,
7

Tagging @sbhasan so that they can see this as well. I believe I have found the other necessary step for interpolating from the parent 2D mesh to 1D submesh which pertains to the parent mesh’s lower boundary. Instead of interpolate, I needed to use interpolate_nonmatching which Dokken provided code for an example regarding a 3D to 2D case HERE. Putting everything together, and also testing the interpolation of Fx onto the original parent mesh’s lower boundary via dirichletbc as this is something that I need to do for my problem, I have the following code tested and running on Dolfinx version 0.9.0

from mpi4py.MPI import COMM_WORLD as comm
import numpy as np
import ufl
from dolfinx.fem import assemble_scalar,dirichletbc,create_interpolation_data,form,Function,functionspace, locate_dofs_geometrical
from dolfinx import plot
from dolfinx.mesh import (
    create_submesh,
    create_unit_square,
    create_rectangle,
    exterior_facet_indices,
    locate_entities,
    locate_entities_boundary,
)
import pyvista as pv
import matplotlib.pyplot as plt
Nx, Ny = 5,5
TYPE = "Lagrange"
ORDER = 2
LABEL_DOFS = "NO"

mesh = create_rectangle(comm,[[0.0,0.0],[np.pi,np.pi]],[Nx,Ny])
tdim = mesh.topology.dim
mesh.topology.create_connectivity(tdim-1,tdim)

def ref_func(x):
    return np.cos(x[0]) + np.sin(x[1]) #This will be cos(x) on the lower boundary


V = functionspace(mesh, (TYPE, ORDER))
f = Function(V)
f.interpolate(ref_func)

# Visualize interpolation of original onto domain

grid = pv.UnstructuredGrid(*plot.vtk_mesh(V))
grid.point_data["f"] = f.x.array.real
grid.set_active_scalars("f")
plotter = pv.Plotter()
plotter.add_mesh(grid, show_edges=True)

if LABEL_DOFS == "YES":
    dof_coords = V.tabulate_dof_coordinates()
    dof_indices = V.dofmap.index_map.local_to_global(np.arange(V.dofmap.index_map.size_local))
    dof_labels = [str(dof_index) for dof_index in dof_indices]

    plotter.add_points(dof_coords,
        color="blue",
        render_points_as_spheres=False,
        point_size=10,
        label="DOFs",
    )
    plotter.add_point_labels(dof_coords,
        dof_labels,
        font_size=10,
        point_color="blue",
        point_size=10,
        always_visible=True,
    )
    
plotter.view_xy()
plotter.show_grid()
plotter.show()
plotter.screenshot("Interpolated_2D_expr.png")

# Locate boundary entities
def lower_boundary_locator(x):
    truth_arr = np.isclose(x[1],0)
    return truth_arr
boundary_cells = locate_entities_boundary(mesh, tdim-1, lower_boundary_locator)

# Create submesh on the boundary
boundary_submesh, parent_cells, vertex_max, node_map = create_submesh(mesh, tdim-1, boundary_cells) # submesh and parent_cells

# Create functionspace using boundary submesh
V_sub_boundary = functionspace(boundary_submesh, (TYPE, ORDER))
fx = Function(V_sub_boundary)

# interpolate reference function onto sub-domain
fine_mesh_cell_map = boundary_submesh.topology.index_map(boundary_submesh.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc,dtype = np.int32)
interpolation_data = create_interpolation_data(V_sub_boundary,V,np.arange(num_cells_on_proc,dtype = np.int32),1e-14)

fx.interpolate_nonmatching(f,cells,interpolation_data = interpolation_data)

# Cumulatively integrate fx

x = ufl.SpatialCoordinate(boundary_submesh)
Fx = Function(V_sub_boundary)
Fx_integral = lambda x0: np.hstack([assemble_scalar(form(fx*ufl.conditional(x[0] < x_, 1, 0)*ufl.dx)) for x_ in x0[0]])
Fx.interpolate(Fx_integral)


# Interpolate Fx onto parent mesh's boundary
boundary_dofs = locate_dofs_geometrical(V,lower_boundary_locator)
bc_lower = dirichletbc(Fx, boundary_dofs)

# Test expression interpolated onto boundary
fx_on_boundary = f.x.array[boundary_dofs]
Fx_interpolated_onto_boundary = bc_lower.g.x.array
x_axis = V.tabulate_dof_coordinates()[boundary_dofs][:,0]
sort_idx = np.argsort(x_axis)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x_axis[sort_idx], fx_on_boundary[sort_idx])
ax.plot(x_axis[sort_idx], Fx.x.array[sort_idx], marker = 'o')
ax.plot(x_axis[sort_idx], Fx_interpolated_onto_boundary[sort_idx])
ax.set_xlabel('$x$')
ax.legend(['$f|_{y=0}$','$Fx_{V_{sub}}$','$Fx_{V_{boundary}}$'])
plt.show()

Looks correct. Thank you all for your help!

1 Like