All of the sudden getting the error AttributeError: 'NoneType' object has no attribute 'ufl_id' when code worked last week

General
1

I was using this code last week that ran and gave me my desired results.

# @title coupler - 3D FEM, 1D dG

"""
The 1D and 3D domains are computed and coupled here.
"""


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_dg(Lambda, file_path, radius_map, uh1d_dg=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_dg(point) for point in points])

    if uh1d_dg 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()

def export_all_for_paraview(Lambda, Omega, uh1d_dg, uh3d_dg, radius_map, output_dir="vtk_output"):
    os.makedirs(output_dir, exist_ok=True)

    # -------------------------- #
    # --- Export 1D Vessel ---  #
    # -------------------------- #
    import pyvista as pv
    coords_1d = Lambda.coordinates()

    poly = pv.PolyData()
    poly.points = coords_1d

    # Build polyline segments from graph edges
    connectivity = []
    for i, j in ves_Geom.edges():
        connectivity.extend([2, i, j])  # VTK format: [n_points, pt0, pt1]

    poly.lines = np.array(connectivity)

    # Point data: pressure and radius
    pressure_1d = np.array([uh1d_dg(p) for p in coords_1d])
    radii_1d = np.array([radius_map(p) for p in coords_1d])

    poly["Pressure1D"] = pressure_1d
    poly["radius"] = radii_1d

    mesh1d_path = os.path.join(output_dir, "lambda_1d.vtp")
    poly.save(mesh1d_path)
    print(f"[EXPORT] Saved 1D vessel to: {mesh1d_path}")

    # -------------------------- #
    # --- Export 3D Tissue ---  #
    # -------------------------- #
    points_3d = Omega.coordinates()
    # Extract tetrahedral connectivity
    tetra_cells = [cell.entities(0) for cell in cells(Omega)]
    cells_3d = [("tetra", np.array(tetra_cells))]

    # Get pressure values at vertices
    pressure_3d = uh3d_dg.compute_vertex_values(Omega)

    mesh3d = meshio.Mesh(
        points_3d,
        cells_3d,
        point_data={"Pressure3D": pressure_3d}
    )
    mesh3d_path = os.path.join(output_dir, "omega_3d.vtu")
    mesh3d.write(mesh3d_path)

    print(f"[EXPORT] Saved 1D mesh to {mesh1d_path}")
    print(f"[EXPORT] Saved 3D mesh to {mesh3d_path}")

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_dg(ves_Geom, directory_path, Lambda, Omega, mf, del_Omega, perf3, perf1, kappa, gamma_in, gamma_out, 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
    """

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

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

    c = Omega.coordinates()
    zmin = np.min(c[:, 2])
    zmax = np.max(c[:, 2])
    zl = zmax - zmin

    #def boundary_Omega(x, on_boundary):
    #    return on_boundary and not near(x[2], zmin) and not near(x[2], zmax)

    # Constants
    kappa = Constant(kappa) #permeability or transfer coefficient
    gamma_in = Constant(gamma_in) #1D boundary coeff
    gamma_out = Constant(gamma_out) #1D boundary coeff
    P = Constant(P) #1D pressure value to be imposed
    #del_Omega = Constant(del_Omega) #3D pressure value to be imposed, comment back in when not doing convergence


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

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

    # Mesh-related terms
    n = FacetNormal(Lambda)
    h = CellDiameter(Lambda)
    alpha = Constant(10.0)  # penalty parameter


    # 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

    #For dsLambda to be used in dG
    boundary_markers = MeshFunction("size_t", Lambda, Lambda.topology().dim() - 1, 0)

    #vessel inlet
    class BottomEnd(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[2], 0.0)

    #vessl outlet
    class TopEnd(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[2], 1.0)  # set to actual domain length in x[0], or use `np.max(Lambda.coordinates()[:, 0])`

    bottom = BottomEnd()
    top = TopEnd()
    bottom.mark(boundary_markers, 1)
    top.mark(boundary_markers, 2)

    dsLambda = Measure("ds", domain=Lambda, subdomain_data=boundary_markers)

    def inlet_3d(x, on_boundary):
        return on_boundary and near(x[2], 0.0)

    def outlet_3d(x, on_boundary):
        return on_boundary and near(x[2], 1.0)

    bc3d = [DirichletBC(V3, del_Omega, inlet_3d),
            DirichletBC(V3, del_Omega, outlet_3d)]

    # 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)

    DG0 = FunctionSpace(Lambda, "DG", 0)
    D_area_expr = Expression("pi*radius*radius", radius=radius_map, pi=np.pi, degree=2)
    D_perimeter_expr = Expression("2*pi*radius", radius=radius_map, pi=np.pi, degree=2)

    D_area = interpolate(D_area_expr, DG0)
    D_perimeter = interpolate(D_perimeter_expr, DG0)

    # 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
    # dG formulation for 1D block
    a11 = (
      perf1*inner(grad(u1), grad(v1))*D_area*dxLambda
      - perf1*dot(avg(grad(u1)), jump(v1, n))*avg(D_area)*dS
      - perf1*dot(avg(grad(v1)), jump(u1, n))*avg(D_area)*dS
      + perf1*alpha('+')/h('+')*dot(jump(u1, n), jump(v1, n))*avg(D_area)*dS
      + kappa*inner(u1, v1)*D_perimeter*dxLambda
      - gamma_in*inner(u1, v1)*D_area*dsLambda(1) + gamma_out*inner(u1, v1)*D_area*dsLambda(2)
  )

    # Right-hand side
    L0 = inner(Constant(0), v3_avg)*dxLambda
    L1 = inner(Constant(0), v1)*D_area*dxLambda - gamma_in*inner(P, v1)*dsLambda(1) + gamma_out*inner(P, v1)*dsLambda(2)

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

    # Define boundary subdomain for 1D
    class BottomBoundary1D(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[2], 0.0)

    bc1d_dg = DirichletBC(V1, Constant(P), BottomBoundary1D(), method="pointwise")

    W_bcs = [
        bc3d,         # 3D
        [bc1d_dg]     # 1D
    ]

    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_dg, uh1d_dg = wh

    # Rename variables
    uh3d_dg.rename("3D Pressure", "3D Pressure Distribution")
    uh1d_dg.rename("1D Pressure", "1D Pressure Distribution")


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

    #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

       # Save for ParaView
    #export_all_for_paraview(Lambda, Omega, uh1d_dg, uh3d_dg, radius_map, output_dir=directory_path)

    return os.path.join(directory_path, "lambda_1d.vtk"), os.path.join(directory_path, "omega_3d.pvd"), uh1d_dg, uh3d_dg, Lambda, Omega

With the call function

output_file_1d, output_file_3d, uh1d_dg, uh3d_dg, Lambda, Omega = steadystate_dg(ves_Geom, OUTPUT_PATH,
    Lambda, Omega, mf,
    del_Omega=800,
    perf3=9.6e-2,
    perf1=1.45e4,
    kappa=3.0,
    gamma_in=1.0,
    gamma_out = 1.0,
    P=934
)

Logging back in this week I am all of the sudden getting the error

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-6-38f996b25037> in <cell line: 0>()
     30 #)
     31 
---> 32 output_file_1d, output_file_3d, uh1d_dg, uh3d_dg, Lambda, Omega = steadystate_dg(ves_Geom, OUTPUT_PATH,
     33     Lambda, Omega, mf,
     34     del_Omega=800,

6 frames
/usr/local/lib/python3.11/dist-packages/dolfin/fem/assembling.py in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
    214                 if hasattr(op, "ufl_function_space"):
    215                     same_mesh = bool(
--> 216                         same_mesh and op.ufl_function_space().ufl_domain().ufl_id() == dolfin_form.mesh().ufl_domain().ufl_id()
    217                     )
    218 

AttributeError: 'NoneType' object has no attribute 'ufl_id'

Was there an update that occurred that changed the compatibility with certain elements of my setup? I am pretty new to fenics and FEM in general so any tips would be appreciated.

2

Hello,

Indeed, we made small changes causing this error yesterday (sorry for inconvenience). It has been fixed today so your code should work again with the latest version.

3

Thank you for your quick reply. I am still getting the error today, I pulled from the following libraries today:

# @title Install required libraries (try-catch safe)
"""
NOTE: You can save the Docker image state after running this block so that you don't have to run it every time you start a new environment.
"""
import os, re

def replace_in_file(file_path):
    with open(file_path, 'r', encoding='utf-8') as file:
        content = file.read()

    # Replace 'ufl' with 'ufl_legacy'
    content = re.sub(r'\bufl\b', 'ufl_legacy', content)

    with open(file_path, 'w', encoding='utf-8') as file:
        file.write(content)

def process_directory(directory):
    for root, _, files in os.walk(directory):
        for file in files:
            if file.endswith('.py'):
                file_path = os.path.join(root, file)
                replace_in_file(file_path)

# ipywidgets
try:
    import ipywidgets
except ImportError:
    !pip install ipywidgets

# dolfin
try:
    import dolfin
except ImportError:
    !wget "https://fem-on-colab.github.io/releases/fenics-install-release-real.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
    import dolfin

# block
try:
    import block
except ImportError:
    !git clone "https://bitbucket.org/fenics-apps/cbc.block/src/master/"
    !pip install master/

# fenics_ii
try:
    import xii
except ImportError:
    !git clone "https://github.com/MiroK/fenics_ii"
    process_directory("fenics_ii/")
    !pip install fenics_ii/

# vtk
try:
    import vtk
except ImportError:
    !pip install vtk

# graphnics
try:
    import graphnics
except ImportError:
    !git clone "https://github.com/IngeborgGjerde/graphnics"
    !pip install graphnics/

# meshio
try:
    import meshio
except ImportError:
    !pip install meshio

# pyvista
try:
    import pyvista
except ImportError:
    !pip install pyvista

I should mention I am running code on Google Colab. Were any of the links to the repositories changed with this fix? Or should I just allow more time for the fix to take place?

4

Hi @al25,
I have triggered a rebuild of the FEniCS package in FEM on Colab, and the fix mentioned above should now be available there.

Keep in mind that if you have a notebook currently running you’ll need to go to “Runtime → Disconnect and delete runtime” to force it to the download the updated version.

2 Likes