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.