To start, I am on the Dolfinx nightly version which is recognized as version DOLFINx version: 0.10.0.0, and I’m on Python version 3.12.3.
Without getting into any of the specifics of the full problem that I’m working on, I need to evaluate the following cumulative integral
\int_0^x v_h|_{y=0} dx
where v_h pertains to the \hat{\textbf{j}}-component of the vector \textbf{V}_h = u_h \hat{\textbf{i}} + v_h \hat{\textbf{j}}
\textbf{V}_h was evaluated in my code via interpolating a ufl expression. Here is my code for creating \textbf{V}_h
from mpi4py import MPI
import pyvista as pv
import dolfinx as df
import ufl
import gmsh
import numpy as np
Nx,Ny = 60,60 # Number of points in x and y
mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, Nx, Ny)
tdim, fdim = mesh.topology.dim, mesh.topology.dim-1
mesh.topology.create_connectivity(fdim, tdim)
# Define a true functionspace for later usage
S = df.fem.functionspace(mesh, ("Lagrange", 1)) # scalar space for grids and projection of scalar quantities
# Define total wind field
V = df.fem.functionspace(mesh, ("Lagrange", 1, (2,))) # 2D Vector function space
wind = df.fem.Function(V)
wind.interpolate(lambda x: (3*(x[0]-x[1]),3*x[0]+x[1]-2))
# Extract vertical relative vorticity from wind
Z_expr = df.fem.Expression(wind[1].dx(0) - wind[0].dx(1),S.element.interpolation_points)
Z = df.fem.Function(S)
Z.interpolate(Z_expr) # vorticity
# Define boundary condition
boundary_facets = df.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = df.fem.locate_dofs_topological(S, fdim, boundary_facets)
bc = df.fem.dirichletbc(df.default_scalar_type(0),boundary_dofs,S)
# Solve
u = ufl.TrialFunction(S)
v = ufl.TestFunction(S)
a = ufl.inner(ufl.grad(u),ufl.grad(v)) * ufl.dx
L = -Z*v*ufl.dx
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a=a,L=L,bcs = [bc], petsc_options_prefix="Streamfunction_problem",petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
psi_I = problem.solve()
# Extract horizontal divergence from wind field (onto scalar space)
d_expr = df.fem.Expression(ufl.div(wind),S.element.interpolation_points)
d = df.fem.Function(S)
d.interpolate(d_expr) # Evaluate div epression
# Solve
u = ufl.TrialFunction(S)
v = ufl.TestFunction(S)
a = ufl.inner(ufl.grad(u),ufl.grad(v)) * ufl.dx
L = -d*v*ufl.dx
problem = LinearProblem(a=a,L=L,bcs = [bc], petsc_options_prefix="Streamfunction_problem",petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
chi_I = problem.solve()
###########################################################################################
# And now we have V_h, and we're ready to begin working on resolving the boundary values ##
###########################################################################################
V_h = df.fem.Function(V) # Define V_h in our 2D vector space, V
V_h_expr = df.fem.Expression((wind - (ufl.grad(chi_I) + ufl.perp(ufl.grad(psi_I)))), V.element.interpolation_points)
V_h.interpolate(V_h_expr) # Evaluate V_h
So it was interpolated onto a 2D functionspace, V, and now I need to cumulatively integrate the v_h component located on the southern boundary of the square i.e. where y=0 and from 0 to x. To do this, I need to extract from V_h all of its v_h component values and find all of them located on the boundary at y=0 from 0\le x\le 1, make sure that these values which I’ve extracted are properly ordered, and then use something like scipy to perform cumulative integration, or maybe use something built into Dolfinx which could do this computation for me. Can this be done?
If the full problem and additional context are necessary, I can provide it, just let me know, but I feel like for what I’m asking, this should be sufficient. In the code, I solved 2 Poisson Equations with homogeneous boundary conditions, and used the 2 solutions to obtain V_h. Now I need to cumulatively integrate a component of V_h.