Hello Susanne,
I had glanced at your paper, very quickly (most of it is way above my head, I am not a specialist in the domain), and today I glanced at OptiCut’s documentation, which is similar in some ways to the information found in the paper. I think I have understood the gist of it, but not the details.
For example, in my case I use polynomials of degree 4 (2 is minimal for me not to get a vanishing quantity that depends on the 2nd order spatial derivatives of “temp” and “volt”, the 2 unknowns of the 2 coupled PDEs in my problem), and from what I understand, the solution is computed not only at the mesh vertices but also in some interior points of the triangles making up the mesh. So I don’t exactly know how CutFEMx is able to retrieve the position of those points, checks whether they are inside or outside the set_level that the user has defined, and so on. But of course, I am confident that it is able to deal with this case without problem.
Unfortunately I do not have a MWE to share. A MWE with all the features that my full code possesses, such as mixed element space, anisotropy in a physical property and a mesh with 2 regions. My code is “big”, spread over several files and uses an object oriented approach, i.e. I created some classes and helper functions to handle my problem. I don’t have much time and I took advantage of LLM AI to create a MWE. I reviewed the code, modified it a bit, executed it and it looks OK to me. It did modify the weak form (you can see this yourself, I gave my real weak form in my previous post), so the equations it solves might not make entirely physical sense. But the code does work and is very similar to my real code, so if we can take it from there, I think I could apply any suggestion to my real code.
The code includes some post processing that can safely be removed, but it’s just some extra lines making sure everything is correct.
Here’s the Python code:
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, log
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar,
dirichletbc, form, locate_dofs_topological, Expression)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import gmshio
from dolfinx import default_scalar_type
from ufl import (Measure, SpatialCoordinate, TestFunction, TrialFunction,
split, derivative, grad, dot, dx, dS, div, FacetNormal)
from basix.ufl import element, mixed_element
# --- 1. Simulation and Material Parameters ---
# Physical constants and simulation settings
T_cold = 300.0
ΔT = 100.0
Îş = 1.8 # Thermal conductivity
# Finite element degree
deg = 4
# Seebeck tensor components for the anisotropic region
s_xx, s_yy = 1e-5, 1e-3
# --- 2. Mesh Loading and Definition of Measures ---
# Load the mesh and associated markers from the .msh file
# This replaces the TheMesh class
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/full_annulus.msh", MPI.COMM_WORLD, gdim=2)
# Define geometric and integration entities
x = SpatialCoordinate(mesh)
n = FacetNormal(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers) # Redundant but kept for consistency
dS = Measure("dS", domain=mesh, subdomain_data=facet_markers)
# --- 3. Definition of Region-Dependent Material Properties ---
# The mesh has two material regions marked with IDs 20 and 21.
# We define material properties (Seebeck tensor, resistivity, conductivity)
# that have different values in each region.
# Create a discontinuous function space (DG0) for region-dependent properties
# Resistivity (rho) and electrical conductivity (sigma)
rho_space = functionspace(mesh, ("DG", 0))
rho = Function(rho_space)
sigma = Function(rho_space)
three_quarters_cells = cell_markers.find(20)
quarter_cells = cell_markers.find(21)
# Assign values for the "three-quarter" material region (ID 20)
rho.x.array[three_quarters_cells] = np.full_like(three_quarters_cells, 1e0, dtype=default_scalar_type)
sigma.x.array[three_quarters_cells] = np.full_like(three_quarters_cells, 1.0 / 1e0, dtype=default_scalar_type)
# Assign values for the "quarter" material region (ID 21)
rho.x.array[quarter_cells] = np.full_like(quarter_cells, 1.0, dtype=default_scalar_type)
sigma.x.array[quarter_cells] = np.full_like(quarter_cells, 1.0 / 1.0, dtype=default_scalar_type)
# Define the anisotropic Seebeck tensor using a DG0 tensor function space
tensor_el = element(family="DG", cell=str(mesh.ufl_cell()), degree=0, shape=(2, 2))
S_space = functionspace(mesh, tensor_el)
Seebeck_tensor = Function(S_space)
# Define the tensor values for each region
def seebeck_three_quarters(x):
# Isotropic with zero Seebeck effect in this region
tensor = np.array([[0, 0], [0, 0]])
num_points = x.shape[1]
return np.tile(tensor.flatten(), (num_points, 1)).T.reshape(4, num_points)
def seebeck_quarter(x):
# Anisotropic Seebeck tensor in this region
tensor = np.array([[s_xx, 0], [0, s_yy]])
num_points = x.shape[1]
return np.tile(tensor.flatten(), (num_points, 1)).T.reshape(4, num_points)
# Interpolate the functions to define the Seebeck tensor across the domain
Seebeck_tensor.interpolate(seebeck_three_quarters, cells0=three_quarters_cells)
Seebeck_tensor.interpolate(seebeck_quarter, cells0=quarter_cells)
# --- 4. Weak Formulation using Mixed Elements ---
# Define a mixed element function space for temperature (temp) and voltage (volt)
el = element("CG", str(mesh.ufl_cell()), deg)
mel = mixed_element([el, el])
ME = functionspace(mesh, mel)
# Define test functions and the function to hold the solution
u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)
# Define the electric current density J
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
# Weak form for the heat equation (F_T)
F_T = dot(Îş * grad(temp), grad(u)) * dx - dot(rho * J_vector, J_vector) * u * dx
F_V = dot(J_vector, grad(v)) * dx
# Combine into the final weak form
weak_form = F_T + F_V
# --- 5. Boundary Conditions ---
# Boundary IDs from the gmsh file
# Physical Curve("r_inner_quarter", 24)
# Physical Curve("r_outer_quarter", 25)
# Physical Curve("r_inner_three_quarter", 26)
# Physical Curve("r_outer_three_quarter", 27)
T_cold_boundary_curve = 24
T_hot_boundary_curve = 25
# Locate degrees of freedom for temperature on the specified boundaries
facets_cold = facet_markers.find(T_cold_boundary_curve)
facets_hot = facet_markers.find(T_hot_boundary_curve)
dofs_cold = locate_dofs_topological(ME.sub(0), mesh.topology.dim - 1, facets_cold)
dofs_hot = locate_dofs_topological(ME.sub(0), mesh.topology.dim - 1, facets_hot)
# Define and apply Dirichlet (fixed value) boundary conditions for temperature
bc_cold = dirichletbc(PETSc.ScalarType(T_cold), dofs_cold, ME.sub(0))
bc_hot = dirichletbc(PETSc.ScalarType(T_cold + ΔT), dofs_hot, ME.sub(0))
bcs = [bc_cold, bc_hot]
# --- 6. Solver Configuration and Solution ---
# The problem is nonlinear, so we define the Jacobian and use a Newton solver
Jac = derivative(weak_form, TempVolt, dTV)
# Set up the nonlinear problem
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J=Jac)
# Configure the Newton solver
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-14
solver.atol = 1e-14
solver.report = True
solver.max_it = 100
# Configure the linear solver (KSP) within the Newton solver
# Use a direct LU solver (MUMPS) for robustness
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
# Solve the system
print("--- Starting nonlinear solver ---")
n_iter, converged = solver.solve(TempVolt)
if converged:
print(f"Solver converged in {n_iter} iterations.")
else:
print("Solver did not converge.")
# The solution is now in TempVolt. We can split it to get temp and volt.
temp_sol, volt_sol = TempVolt.sub(0).collapse(), TempVolt.sub(1).collapse()
# --- 7. Post-Processing (Optional) ---
# Example calculation: Check for current conservation across internal boundaries
# Physical Curve("theta_0", 22)
# Physical Curve("theta_pi_over_2", 23)
# These are internal lines separating the two material regions
J_integrand = (-sigma('+') * grad(volt)('+') - sigma('+') * Seebeck_tensor('+') * grad(temp)('+'))
# Calculate the total current flowing through each interface
I_theta0 = assemble_scalar(form(dot(J_integrand, n('+')) * dS(22)))
I_thetapi_over_2 = assemble_scalar(form(dot(J_integrand, n('+')) * dS(23)))
print("\n--- Post-processing Results ---")
print(f"Current through interface at theta=0 (ID 22): {I_theta0:.6e} A")
print(f"Current through interface at theta=pi/2 (ID 23): {I_thetapi_over_2:.6e} A")
# For a closed system with no current sources/sinks, these should sum to zero.
# Check that the divergence of J is close to zero over the whole domain
div_J_integral = assemble_scalar(form(div(J_vector) * dx))
print(f"Volume integral of div(J): {div_J_integral:.6e}")
Its execution yields:
Info : Reading 'meshes/full_annulus.msh'...
Info : 13 entities
Info : 6205 nodes
Info : 12444 elements
Info : Done reading 'meshes/full_annulus.msh'
--- Starting nonlinear solver ---
Solver converged in 2 iterations.
--- Post-processing Results ---
Current through interface at theta=0 (ID 22): 7.841666e-03 A
Current through interface at theta=pi/2 (ID 23): 8.043990e-03 A
Volume integral of div(J): -3.967311e-12
The mesh I am considering comes from the following file.geo:
// Gmsh project created on Sat Mar 8 22:37:56 2025
SetFactory("OpenCASCADE");
//+
Mesh.CharacteristicLengthMax = 0.0009;
Mesh.CharacteristicLengthMin = 0.0001;
Disk(1) = {0, 0, 0, 0.05, 0.05};
//+
Disk(2) = {0, 0, 0, 0.035, 0.035};
BooleanDifference{ Surface{1}; Delete;}{ Surface{2}; Delete;}
// Surface 3 is the full annulus. I need to split it into 2.
Point(4) = {0, 0.05, 0};
//+
Point(5) = {0, 0.035, 0};
//+
// Origin.
Point(6) = {0, 0, 0};
Line(4) = {2,3};
Line(5) = {4,5};
Circle(6) = {5, 6, 2};
Circle(7) = {4, 6, 3};
Curve Loop(8) = {5,6,4,7};
Plane Surface(9) = {8};
// 3-quarter surface
BooleanDifference{ Surface{1}; Delete;}{ Surface{9};}
BooleanFragments{ Surface{1}; Surface{9}; Delete; }{ }
Physical Surface("three_quarter_material", 20) = {1};
Physical Surface("quarter_material", 21) = {9};
Physical Curve("theta_0", 22) = {4};
Physical Curve("theta_pi_over_2", 23) = {5};
Physical Curve("r_inner_quarter", 24) = {6};
Physical Curve("r_outer_quarter", 25) = {7};
Physical Curve("r_inner_three_quarter", 26) = {8};
Physical Curve("r_outer_three_quarter", 27) = {9};
which can generate the file.msh either from Gmsh GUI or from the command gmsh -2 full_annulus.geo.
Let me know if that’s too much, too big, i.e. if you would prefer something trimmed. I just don’t know if I would be able to apply the suggestions to my real code if I trim it significantly from the above example.
Thank you very much for your time and your help. 
