Use of cutfemx for a thermoelectric problem

[mantledparamountic]

Dear Susanne and other Cutfemx developper(s), I would really appreciate some help to solve a thermoelectric problem, which I am already able to solve with standard dolfinx 0.9. My end goal would be to modify the level set of a subdomain, to modify a part of the geometry of the material where the PDEs are solved, and perform a shape optimization under some constraints. I already posted some attempts here, but I am facing several problems.

The domain is an annulus, with the first quadrant (quarter-annulus) made of a particular thermoelectric material with anisotropic Seebeck "coefficient", or tensor (a 2x2 diagonal matrix in my case, with different entries). The remaining 3 quarters are made up of an electrically and thermally conducting material but with vanishing Seebeck coefficient. I am solving div J_U =0 and div J = 0 where J = -sigma grad(volt) - sigma S grad(temp), J_U = -kappa grad (temp) + T SJ.
As boundary conditions, I apply Dirichlet b.c.s for the temperature at the inner and outer radii of the annulus, or quarter-annulus or even at the theta=0 and theta=pi/2 interfaces, which are the 2 interfaces between the materials. Each one of those cases represent a different problem I am already able to solve using standard Fenicsx. The solver is able to deal with the voltage's lack of b.c.s. and automatically fixes the 0 of the electrostatics potential (my volt variable).

If we solve the problem using CutFEMx, I would not mind if you took it as an example in the doc or tutorial, it uses mixed elements with 2 different materials/regions and so is a bit more involved than the usual simple examples. That's also a reason why I post my problem here, maybe it is better suited than in the FEniCSx discourse forums.

The very first question I have is whether my level set is correct, and if so, why using a high density mesh (square of 300x300) creates the error

ValueError: triangle is not intersected and therefore cannot be cut

This is counter-intuitive to me. If I relax the mesh density to a square of 200x200 or 50x50, the error goes out.

My code for the level set is here:

import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
import ufl
from ufl import (FacetNormal, CellDiameter,
                 conditional, gt, sqrt,
                 grad, inner, dot, jump, avg)

from dolfinx import fem, mesh
from dolfinx import default_scalar_type
from dolfinx.fem import (functionspace, Function, Constant)


from cutfemx.level_set import (
    locate_entities as cutfemx_locate_entities,
    compute_normal,
)
from cutfemx.quadrature import runtime_quadrature
from cutfemx.petsc import assemble_vector
from cutfemx.fem import cut_form
from basix.ufl import element, mixed_element
from dolfinx.mesh import GhostMode
from dolfinx import plot
import pyvista as pv


T_cold = 300.0
T_hot = 400.0
kappa = 1.8
s_xx, s_yy = 1e-5, 1e-3
deg = 2
gamma_N = 100 * deg**2
R_inner = 0.035
R_outer = 0.05
N = 200

msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=((-0.15, -0.15), (0.15, 0.15)),
    n=(N, N),
    cell_type=mesh.CellType.triangle,
    ghost_mode=GhostMode.shared_facet,
)
tdim = msh.topology.dim # 2.
gdim = msh.geometry.dim # 2.

V_level_set = fem.functionspace(msh, ("Lagrange", deg))

def annulus_level_set_func(x):
    r_val = np.sqrt(x[0]**2 + x[1]**2)
    return np.maximum(r_val - R_outer, R_inner - r_val)

level_set = fem.Function(V_level_set)
level_set.interpolate(annulus_level_set_func)

# Mixed FE space for (temperature, voltage)
el0 = element("CG", msh.ufl_cell().cellname(), deg)
ME  = fem.functionspace(msh, mixed_element([el0, el0]))
UV  = Function(ME)
dUV = ufl.TrialFunction(ME)
v_T, v_V = ufl.TestFunctions(ME)
temp, volt = ufl.split(UV)

# Material + Seebeck in Q1 only
sigma = Constant(msh, default_scalar_type(1.0))
x = ufl.SpatialCoordinate(msh)
r = sqrt(x[0] ** 2 + x[1] ** 2)

in_ann  = ufl.And(r >= R_inner, r <= R_outer)
in_q1   = ufl.And(x[0] >= 0, x[1] >= 0)
marker  = conditional(ufl.And(in_ann, in_q1), 1.0, 0.0)
Smat    = ufl.as_matrix([[s_xx, 0], [0, s_yy]])
Seebeck = marker * Smat



# Plot the mesh, level set and region where S is different from 0.###########################

# only on rank 0
if MPI.COMM_WORLD.rank == 0:
    # 1) Build a VTK‐compatible unstructured grid of V_level_set
    cells, types, coords = plot.vtk_mesh(V_level_set)
    grid = pv.UnstructuredGrid(cells, types, coords)

    # 2) Attach the level‐set φ as a point‐scalar
    grid.point_data["phi"] = level_set.x.array

    # 3) (optional) warp by φ to see the "height" of φ
    warped = grid.warp_by_scalar(scalars="phi", factor=0.01)

    # 4) Plot the colored field, the warped mesh, and φ=0 contour
    p = pv.Plotter(window_size=[800, 600])
    p.add_mesh(
        grid,
        scalars="phi",
        cmap="coolwarm",
        show_scalar_bar=True,
        lighting=False,
        opacity=0.6,
        name="level_set"
    )
    p.add_mesh(
        warped,
        show_edges=True,
        color="white",
        opacity=0.4,
        name="warp"
    )
    # zero‐level contour on the original grid
    contour = grid.contour([0.0], scalars="phi")
    p.add_mesh(
        contour,
        color="black",
        line_width=3,
        name="phi=0"
    )

    p.camera_position = "xy"    # 2D top‐down view
    p.add_text("Level‐set φ(x,y)", font_size=14, color="black")
    p.show()


# extract negative and positive regions
neg = grid.threshold(0.0, scalars="phi", invert=True)   # φ < 0
pos = grid.threshold(0.0, scalars="phi", invert=False)  # φ ≥ 0

p = pv.Plotter()

# 1) φ < 0 → solid red
p.add_mesh(neg, color="red", show_scalar_bar=False)

# 2) φ ≥ 0 → solid blue
p.add_mesh(pos, color="blue", show_scalar_bar=False)

# 3) zero‐level contour → white (or black)
contour = grid.contour([0.0], scalars="phi")
p.add_mesh(contour, color="white", line_width=4)

p.add_mesh(grid, show_edges=True, edge_color="black", opacity=0.3)


p.camera_position = "xy"
p.add_text("φ < 0 (red) | φ > 0 (blue) | φ=0 (white)", font_size=12)
p.show()



# 1) Build a CG1 space
V_S = functionspace(msh, ("Lagrange", 1))
S_xx = Function(V_S)

# 2) Analytical Seebeck‐marker on (annulus ∩ Q1)
def seebeck_marker(x):
    # x is shape (2, Npoints)
    r = np.sqrt(x[0]**2 + x[1]**2)
    # φ < 0 ⇔ inside the annulus
    phi = np.maximum(r - R_outer, R_inner - r)
    # mask = inside annulus AND in Q1
    mask = np.logical_and.reduce([phi < 0,
                                  x[0] >= 0,
                                  x[1] >= 0])
    # only S_xx on that region
    return mask.astype(np.float64) * s_xx

S_xx.interpolate(seebeck_marker)

# 3) Extract VTK mesh + coords
cells, types, coords = plot.vtk_mesh(V_S)

# 4) Create PyVista grid and attach the S_xx data
sgrid = pv.UnstructuredGrid(cells, types, coords)
sgrid.point_data["S_xx"] = S_xx.x.array

# 5) Plot it
p = pv.Plotter()
p.add_mesh(
    sgrid,
    scalars="S_xx",
    cmap="hot",
    show_edges=True,
    show_scalar_bar=True,
    lighting=False,
)
p.add_text("Seebeck $S_{xx}$ on (annulus ∩ Q1)", font_size=14)
p.camera_position = "xy"
p.show()
##################################################################################

Whereas my full attempt at solving the problem is here:

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from ufl import (
    FacetNormal, CellDiameter, grad, inner, dot, jump, avg,
    split, derivative, TrialFunction, TestFunctions, conditional, And, sqrt,
    as_matrix, SpatialCoordinate
)

from dolfinx import fem, mesh, log
from dolfinx.fem import (
    functionspace, Function, dirichletbc
)
# FINAL FIX: Explicitly import the petsc submodule to access its functions
import dolfinx.fem.petsc
from dolfinx.mesh import GhostMode
from basix.ufl import element, mixed_element

from cutfemx.level_set import ghost_penalty_facets
from cutfemx.quadrature import runtime_quadrature
from cutfemx.petsc import assemble_vector, assemble_matrix, locate_dofs
from cutfemx.fem import cut_form


# --- 1. Simulation and Material Parameters ---
T_cold = 300.0
T_hot = 400.0
kappa = 1.8
deg = 2
s_xx, s_yy = 1e-5, 1e-3
gamma_N = 10 * deg**2
gamma_g = 1e-3
R_inner = 0.035
R_outer = 0.05
points = ((-0.15, -0.15), (0.15, 0.15))
N = 50
max_iter = 25
rtol = 1e-9

# --- 2. Background Mesh and Level Set Definition ---
msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=points,
    n=(N, N),
    cell_type=mesh.CellType.triangle,
    ghost_mode=GhostMode.shared_facet,
)

V_level_set = functionspace(msh, ("Lagrange", deg))
def annulus_level_set_func(x):
    r_val = np.sqrt(x[0]**2 + x[1]**2)
    return np.maximum(r_val - R_outer, R_inner - r_val)

level_set = Function(V_level_set)
level_set.interpolate(annulus_level_set_func)

# --- 3. Definition of Measures and Function Spaces ---
order = 2 * deg
q_inside = runtime_quadrature(level_set=level_set, ls_part="phi<0", order=order)
q_interface = runtime_quadrature(level_set=level_set, ls_part="phi=0", order=order)

dx_c = ufl.Measure("dC", subdomain_data=[(0, q_inside)], domain=msh)
dC = ufl.Measure("dC", subdomain_data=[(1, q_interface)], domain=msh)

gp_facets = ghost_penalty_facets(level_set, "phi<0")
dS_g = ufl.Measure("dS", subdomain_data=[(0, gp_facets)], domain=msh)

el0 = element("CG", msh.ufl_cell().cellname(), deg)
ME = functionspace(msh, mixed_element([el0, el0]))

# --- 4. Stabilizing Dirichlet BC for Exterior Domain ---
V0, _ = ME.sub(0).collapse()
V1, _ = ME.sub(1).collapse()
exterior_dofs_T = locate_dofs("phi>0", level_set, [V0])
exterior_dofs_V = locate_dofs("phi>0", level_set, [V1])

bc_exterior_T = dirichletbc(PETSc.ScalarType(0), exterior_dofs_T, V0)
bc_exterior_V = dirichletbc(PETSc.ScalarType(0), exterior_dofs_V, V1)
bcs = [bc_exterior_T, bc_exterior_V]

# --- 5. Weak Formulation using CutFEM ---
TempVolt = Function(ME)
(temp, volt) = split(TempVolt)
(v_T, v_V) = TestFunctions(ME)
dTV = TrialFunction(ME)

x = SpatialCoordinate(msh)
r = sqrt(x[0]**2 + x[1]**2)
h = CellDiameter(msh)

in_q1 = And(x[0] >= 0, x[1] >= 0)
marker_q1 = conditional(in_q1, 1.0, 0.0)

rho = marker_q1 * 1.0 + (1 - marker_q1) * 1e0
sigma = 1.0 / rho
Smat = as_matrix([[s_xx, 0], [0, s_yy]])
Seebeck_tensor = marker_q1 * Smat

J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)

F_T = dot(kappa * grad(temp), grad(v_T)) * dx_c
F_T += dot(temp * Seebeck_tensor * J_vector, grad(v_T)) * dx_c
F_T += dot(volt * J_vector, grad(v_T)) * dx_c
F_V = -dot(grad(v_V), sigma * grad(volt)) * dx_c
F_V -= dot(grad(v_V), sigma * Seebeck_tensor * grad(temp)) * dx_c

n_ls = grad(level_set) / sqrt(dot(grad(level_set), grad(level_set)) + 1e-8)
T_g = conditional(r < (R_inner + R_outer) / 2, T_cold, T_hot)
marker_bc_q1 = conditional(in_q1, 1.0, 0.0)
nitsche_term = - dot(kappa * grad(temp), n_ls) * v_T
nitsche_term += - dot(kappa * grad(v_T), n_ls) * (temp - T_g)
nitsche_term += (gamma_N / h) * kappa * (temp - T_g) * v_T
F_T += marker_bc_q1 * nitsche_term * dC

n = FacetNormal(msh)
ghost_penalty = avg(gamma_g) * avg(h) * dot(jump(grad(temp), n), jump(grad(v_T), n)) * dS_g(0)
ghost_penalty += avg(gamma_g) * avg(h) * dot(jump(grad(volt), n), jump(grad(v_V), n)) * dS_g(0)

weak_form_ufl = F_T + F_V + ghost_penalty
Jac_ufl = derivative(weak_form_ufl, TempVolt, dTV)

# --- 6. Create cut_forms ---
F_cut = cut_form(weak_form_ufl)
J_cut = cut_form(Jac_ufl)

# --- 7. Solver Configuration ---
opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_mat_solver_type"] = "mumps"

solver = PETSc.KSP().create(msh.comm)
solver.setFromOptions()

dTV_sol = Function(ME)
TempVolt.x.array[:] = 0.0

# --- 8. Manual Newton Solver ---
log.set_log_level(log.LogLevel.INFO)
print("--- Starting nonlinear solver ---")

for k in range(max_iter):
    A = assemble_matrix(J_cut, bcs=bcs)
    A.assemble()

    b = assemble_vector(F_cut)
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    
    # Apply BCs to residual vector using the fully qualified name
    dolfinx.fem.petsc.set_bc(b, bcs)
    b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

    residual_norm = b.norm()
    if k == 0:
        initial_residual_norm = b.norm()
        if np.isnan(initial_residual_norm):
            print("Solver failed on first iteration: residual is nan.")
            break
    
    print(f"Iteration {k:2d}: Residual norm = {residual_norm:.3e}")
    if residual_norm < rtol * initial_residual_norm:
        print("Solver converged.")
        break

    solver.setOperators(A)
    solver.solve(-b, dTV_sol.vector)
    
    dTV_sol.x.scatter_forward()
    TempVolt.x.array[:] += dTV_sol.x.array
    TempVolt.x.scatter_forward()

if k == max_iter - 1:
    print("Solver did not converge.")

# --- 9. Post-Processing and Visualization ---
from dolfinx.io import VTXWriter
temp_sol, volt_sol = TempVolt.sub(0).collapse(), TempVolt.sub(1).collapse()
temp_sol.name = "Temperature"
volt_sol.name = "Voltage"

with VTXWriter(msh.comm, "annulus_thermoelectric_final.bp", [temp_sol, volt_sol]) as vtx:
    vtx.write(0.0)

print("\nSolution saved to annulus_thermoelectric_final.bp")

I think I get nan's due to my ghost term penalty, but I do not know what's wrong with it.

[mantledparamountic]

If I add a plot to display quadrature points, based on the demo example, I get:

Would this seem "correct"? I see quadrature points that are inside triangles that are not directly on the boundaries of the annulus, though they are the 1st neighbors of those on the boundaries. I suppose this is fine, but I don't know enough about finite elements to be sure. I would like to be sure everything's fine.

The code is

# =================================================================
# Visualize Quadrature Points inside the Annulus (phi < 0)
# This section replicates the behavior of demo_runtime_quadrature.py
# to display the quadrature points on a mesh fitted to the annulus.
# =================================================================

# The region of interest for volume integration is the annulus, where level_set < 0
ls_tag = "phi<0"

# 1. Locate cells fully inside the annulus and cells intersected by the boundary
inside_cells = cutfemx_locate_entities(level_set, tdim, ls_tag)
intersected_cells = cutfemx_locate_entities(level_set, tdim, "phi=0")

# 2. Create a new mesh that is fitted to the annulus geometry.
# This involves creating new sub-cells for the elements cut by the level set.
dof_coords = V_level_set.tabulate_dof_coordinates()
cut_cells_info = cut_entities(level_set, dof_coords, intersected_cells, tdim, ls_tag)
# The cut mesh is composed of the original cells fully inside the annulus
# and the newly created cells that are cut by the boundary.
cut_msh = create_cut_mesh(msh.comm, cut_cells_info, msh, inside_cells)


# 3. Define a runtime quadrature rule for the volume of the annulus
quad_order = 2 #* deg
q_vol = runtime_quadrature(level_set=level_set, ls_part=ls_tag, order=quad_order)

# 4. Compute the physical locations of the quadrature points.
# The `q_vol` object contains the information about which cells to integrate over.
# `physical_points` maps the reference quadrature points to the physical elements.
points_phys = physical_points(q_vol, msh)


# 5. Plot the cut mesh and the quadrature points using PyVista
if MPI.COMM_WORLD.rank == 0:
    plotter = pv.Plotter(window_size=[800, 800])

    # Create a PyVista grid from the newly created cut_mesh
    cells, types, x = plot.vtk_mesh(cut_msh._mesh)
    grid = pv.UnstructuredGrid(cells, types, x)

    # Add the cut mesh to the plotter
    plotter.add_mesh(grid, show_edges=True, color='lightgray', opacity=0.7,
                     label="Fitted Mesh")

    # Add the quadrature points to the plotter
    plotter.add_points(points_phys, render_points_as_spheres=True,
                       point_size=5, color='darkred', label="Quadrature Points")

    # Add the zero-level contour of the original geometry for reference
    cells_orig, types_orig, coords_orig = plot.vtk_mesh(V_level_set)
    original_grid = pv.UnstructuredGrid(cells_orig, types_orig, coords_orig)
    original_grid.point_data["phi"] = level_set.x.array
    contour = original_grid.contour([0.0], scalars="phi")
    plotter.add_mesh(contour, color="black", line_width=3, label="Annulus Boundary")

    plotter.view_xy()
    plotter.add_text(f"Quadrature Points (order {quad_order}) in Annulus", font_size=14)
    plotter.add_legend()
    plotter.show()

which I added directly below/after the code I posted above, related to the level_set.

[mantledparamountic]

Another question I have is whether my problem is unsolvable using CutFEMx due to the use of mixed elements? Is there a remedy I could use if this was the case?

[mantledparamountic]

Hello Susanne (and other devs!), I have seen that you created a branch for a compatibility with mixed elements. This is obviously great, but I wanted to let you know I could "bypass" my problem using standard FEniCSx, with a standard approach of mesh optimization. I.e. I had to remesh a new domain in every cases I considered in my optimization. It was a bit costly but not prohibitive, so I could actually "finish" my work and answer the questions I had.

So, no hurry at all from my part. I might give cutfemx a try in the future, for sure, but I don't "absolutely" need it as soon as possible. Just letting you know.