Help computing exact solution

Hi I am new to Fenics and also new to FEM. Here I am trying to test my code by verifying with an exact solution:

import vtk
import meshio

from dolfin import *
from vtk.util.numpy_support import vtk_to_numpy
from xii import *
from graphnics import *
from scipy.spatial import cKDTree

def save_mesh_vtk(Lambda, file_path, radius_map, uh1d=None):
    """
    Saves a tube mesh as a VTK file with the option to include the 1D solution as a data array.

    Args:
        Lambda (dolfin.Mesh): The mesh to be saved.
        file_path (str): The path where the VTK file will be saved.
        radius_map (dolfin.UserExpression): A function to compute radius values at each point.
        uh1d (dolfin.Function, optional): A function representing 1D pressure data.

    Returns:
        None
    """
    points = Lambda.coordinates()
    cells = {"line": Lambda.cells()}

    radius_values = np.array([radius_map(point) for point in points])
    uh1d_values = np.array([uh1d(point) for point in points])

    if uh1d is not None:
        mesh = meshio.Mesh(points, cells, point_data={"radius": radius_values, "Pressure1D": uh1d_values})
    else:
        mesh = meshio.Mesh(points, cells, point_data={"radius": radius_values})
    mesh.write(file_path)

    reader = vtk.vtkUnstructuredGridReader()
    reader.SetFileName(file_path)
    reader.Update()

    geometryFilter = vtk.vtkGeometryFilter()
    geometryFilter.SetInputData(reader.GetOutput())
    geometryFilter.Update()

    polydata = geometryFilter.GetOutput()

    writer = vtk.vtkPolyDataWriter()
    writer.SetFileName(file_path)
    writer.SetInputData(polydata)
    writer.Write()

class radius_function(UserExpression):
    """
    A user expression to compute the radius at a given point based on the nearest control point in a graph.

    Args:
        ves_Geom (graphnics.FenicsGraph): The graph containing control points with radius information.
        mf (dolfin.MeshFunction): A mesh function associated with the graph.
        kdtree (scipy.spatial.cKDTree): A k-d tree for efficient nearest neighbor search in the graph.
        **kwargs: Additional keyword arguments to be passed to the UserExpression constructor.
    """
    def __init__(self, ves_Geom, mf, kdtree, **kwargs):
        self.ves_Geom = ves_Geom
        self.mf = mf
        self.kdtree = kdtree
        super().__init__(**kwargs)

    def eval(self, value, x):
        p = (x[0], x[1], x[2])
        _, nearest_control_point_index = self.kdtree.query(p)
        nearest_control_point = list(self.ves_Geom.nodes)[nearest_control_point_index]
        value[0] = self.ves_Geom.nodes[nearest_control_point]['radius']

    def value_shape(self):
        return ()

def steadystate(ves_geom, directory_path, del_Omega, perf3, perf1, kappa, gamma, P):
    """
    Runs a steady state simulation with prescribed pressure in both the 1D and 3D domains

    Args:
        G (graphnics.FenicsGraph): The graph representing the network.
        directory_path (str): The directory where the results will be saved
        del_Omega (float, optional): Boundary condition value for 3D pressure. Defaults to 3.0
        perf3     (float, optional): 3D perfusion coefficient
        perf1     (float, optional): 1D perfusion coefficient
        kappa     (float, optional): Transfer coefficient, defines the permeability
        gamma     (float, optional): Boundary condition coefficient for 1D pressure
        P         (float, optional): Boundary condition value for 1D pressure

    Returns:
        tuple (an ordered, unchanging list): A tuple containing:
            - output_file_1d (str): The path to the saved 1D pressure VTK file
            - output_file_3d (str): The path to the saved 3D pressure PVD file
            - uh1d (dolfin.Function): The computed 1D pressure function
            - uh3d (dolfin.Function): The computed 3D pressure function
    """
    # Create \Lambda
    ves_Geom.make_mesh()
    Lambda, mf = ves_Geom.get_mesh()

    # Reference copy of \Lambda
    H = ves_Geom.copy()

    # Create \Omega
    Omega = UnitCubeMesh(20, 10, 10)

    pos = nx.get_node_attributes(ves_Geom, "pos")
    node_coords = np.asarray(list(pos.values()))
    xmin, ymin, zmin = np.min(node_coords, axis=0)
    d = Lambda.coordinates()
    d[:, :] += [-xmin, -ymin, -zmin]
    for node in H.nodes:
        H.nodes[node]['pos'] = np.array(H.nodes[node]['pos']) + [-xmin, -ymin, -zmin]

    kdtree = cKDTree(np.array(list(nx.get_node_attributes(H, 'pos').values())))

    c = Omega.coordinates()
    xl, yl, zl = (np.max(node_coords, axis=0)-np.min(node_coords, axis=0))
    c[:,:] *= [xl+3, yl+3, zl]

    def boundary_Omega(x, on_boundary):
        return on_boundary and not near(x[2], 0) and not near(x[2], zl)

    # Constants
    kappa = Constant(kappa) #permeability or transfer coefficient
    gamma = Constant(gamma) #1D boundary coeff
    P = Constant(P) #1D pressure value to be imposed
    del_Omega = Constant(del_Omega) #3D pressure value to be imposed

    # Function spaces
    V1 = FunctionSpace(Lambda, "CG", 1)
    V3 = FunctionSpace(Omega, "CG", 1)
    W = [V3, V1]

    u3, u1 = list(map(TrialFunction, W))
    v3, v1 = list(map(TestFunction, W))


    # define radius values for vessel - modeled as cylinder
    radius_map = radius_function(ves_Geom, mf, kdtree)
    cylinder = Circle(radius=radius_map, degree=5)

    # Define average of the cross-section
    u3_avg = Average(u3, Lambda, cylinder)
    v3_avg = Average(v3, Lambda, cylinder)

    # Dirac measures
    dxOmega = Measure("dx", domain=Omega) #3D
    dxLambda = Measure("dx", domain=Lambda) #1D
    dsLambda = Measure("ds", domain=Lambda) #1D boundary

    # Define D_area and D_perimeter
    D_area = Expression("pi*radius*radius", radius=radius_map, pi=np.pi, degree=2) # parameterization of cross-section
    D_perimeter = Expression("2*pi*radius", radius=radius_map, pi=np.pi, degree=2) # boundary of cross-section (circumference)



    # Blocks
    a00 = perf3*inner(grad(u3), grad(v3))*dxOmega + kappa*inner(u3_avg, v3_avg)*D_perimeter*dxLambda
    a01 = -kappa*inner(u1, v3_avg)*D_perimeter*dxLambda
    a10 = -kappa*inner(u3_avg, v1)*D_perimeter*dxLambda
    a11 = perf1*inner(grad(u1), grad(v1))*D_area*dxLambda + kappa*inner(u1, v1)*D_perimeter*dxLambda - gamma*inner(u1, v1)*dsLambda

    # Right-hand side
    L0 = inner(Constant(0), v3_avg)*dxLambda
    L1 = inner(Constant(0), v1)*D_area*dxLambda - gamma*inner(P, v1)*dsLambda

    a = [[a00, a01], [a10, a11]]
    L = [L0, L1]

    W_bcs = [[DirichletBC(V3, del_Omega, boundary_Omega)], []]

    A, b = map(ii_assemble, (a, L))
    A, b = apply_bc(A, b, W_bcs)
    A, b = map(ii_convert, (A, b))

    wh = ii_Function(W)

    solver = LUSolver(A, "mumps")
    solver.solve(wh.vector(), b)
    uh3d, uh1d = wh

    # Rename variables
    uh3d.rename("3D Pressure", "3D Pressure Distribution")
    uh1d.rename("1D Pressure", "1D Pressure Distribution")

    # --------------------- #
    # VERIFY EXACT SOLUTION #
    # --------------------- #

    u3e = ii_Function(V3)
    u3e.assign(Constant(P_val))

    u1e = ii_Function(V1)
    u1e.assign(Constant(P_val))

    # Now u3e carries mesh info → Average will work
    u3e_avg = Average(u3e, Lambda, cylinder)

    # Test functions
    v3 = TestFunction(V3)
    v1 = TestFunction(V1)

    # u3e is constant, but we still need to lift it to 1D domain for coupling
    u3e_avg = Average(u3e, Lambda, cylinder)

    # Residuals using constants
    residual_3d = ii_assemble(perf3*inner(grad(u3e), grad(v3))*dxOmega +
                              kappa*inner(u3e_avg, v3)*D_perimeter*dxLambda -
                              kappa*inner(u1e, v3)*D_perimeter*dxLambda)

    residual_1d = ii_assemble(perf1*inner(grad(u1e), grad(v1))*D_area*dxLambda +
                              kappa*inner(u1e, v1)*D_perimeter*dxLambda -
                              kappa*inner(u3e_avg, v1)*D_perimeter*dxLambda -
                              gamma*inner(u1e, v1)*dsLambda -
                              gamma*inner(P, v1)*dsLambda)

    print(f"Residual norm (3D): {residual_3d.norm('l2')}")
    print(f"Residual norm (1D): {residual_1d.norm('l2')}")


    # ------------------------- #
    # SAVE AND RETURN SOLUTIONS #
    # ------------------------- #

    # Create output directory if it doesn't exist and save
    os.makedirs(directory_path, exist_ok=True)
    output_file_1d = os.path.join(directory_path, "pressure1d.vtk")
    output_file_3d = os.path.join(directory_path, "pressure3d.pvd")
    save_mesh_vtk(Lambda, output_file_1d, radius_map, uh1d=uh1d)
    File(output_file_3d) << uh3d

    return output_file_1d, output_file_3d, uh1d, uh3d, Lambda, Omega

But I am getting the error:

TypeError Traceback (most recent call last)
in <cell line: 0>()
4 # pascal transformation from mmHg for pressures
5 # This will store uh1d, uh3d, Lambda, and Omega in RAM
----> 6 output_file_1d, output_file_3d, uh1d, uh3d, Lambda, Omega = steadystate(ves_Geom, OUTPUT_PATH,
7 del_Omega=3.0,
8 perf3=9.6e-2,

1 frames
/usr/local/lib/python3.11/dist-packages/xii/linalg/function.py in init(self, W, components)
21 def init(self, W, components=None):
22 if components is None:
—> 23 self.functions = list(map(Function, W))
24 else:
25 assert len(components) == len(W)

TypeError: ‘FunctionSpace’ object is not iterable

If I remove the section where I am verifying the solution, the code works fine. But I am very interested in verifying the solution and looking for a way to make it work. It seems when I debug one issue then I run into another error so maybe I am not familiar enough with the innerworkings of fenics.

I should also mention, I am running this to test on a straight line:

x1, y1, z1 = (0, 0, 0)
x2, y2, z2 = (2, 0, 2)

ves_Geom = FenicsGraph()
ves_Geom.add_node(0, pos=(x1, y1, z1), radius = 0.5, damage = 1.0)
ves_Geom.add_node(1, pos=(x2, y2, z2), radius = 0.5, damage = 1.0)

ves_Geom.add_edge(0, 1)
print(ves_Geom)