How can I implement "complete electrode model" in dolfinx

variational formulation
1

I am brand new to fenics and I am trying to implement the following strong form:
image
zsh_\ell, sigma, and i_\ell are all known.

I have this as a weak from:
image

In particular, there are two extra (not in the mesh) degrees of freedom here, v_1 and v_2. I discretized it as follows:
image
image

Where the bottom block row has the two extra equations for the two extra dofs.

This is based on the complete electrode model (see for example https://doi.org/10.1137/015206 or https://doi.org/10.3934/math.2021431)

Note that for my application, sigma is a 3x3 tensor, but I think i figured out how to that in dolfinx by just declaring, for example,

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)
basis = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2))
u = ufl.TrialFunction(basis)
v = ufl.TestFunction(basis)
sigma = dolfinx.fem.Constant(
    mesh,
    dolfinx.default_scalar_type([
        [1,0,0],
        [0,1,0],
        [0,0,1],
]))
A = ufl.inner(sigma * ufl.grad(u), ufl.grad(v)) * ufl.dx
# Don't know how to get B,C,D above
M = # Don't know how to get full system matrix (block matrix?)
a = # or maybe I'm trying to get just a lhs representation?
f = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1))
L = ufl.inner(f, v) * ufl.dx
...
bc = # don't know what to put here
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

Following the tutorials, it was pretty straight forward to get the above code. But I cannot figure out how to assemble the other terms over only the respective e_\ell domains or how to combine all the terms into one LinearProblem to solve. I’ve written it as a block matrix, but I’m not familiar enough with doxfinx terms or running on a cluster (which I am) to know if this is best approached as a block matrix type of problem. Can you point me to some examples of this kind of assembly?

I am targeting an HPE Cray machine where I have dolfinx installed (dolfinx.__version__=='.0.7.2') and working via micromamba and conda-forge. My problem scope is 100+ time-harmonic frequencies, each an independent problem with 10+ million dofs. I mention this in case there are pitfalls waiting for me when I try to scale up from toy problems.

2

I think I’m making some progress here… this tutorial has some example of how to integrate over subdomains: Hyperelasticity — FEniCSx tutorial

My domain looks like a solid box with two cylinders cut from the interior. The exterior surface of the box gets a dirichlet condition at zero. The surface of each interior hole is an “electrode” or e_\ell above.

So I am locating the electrodes and wall geometrically using this tutorial Implementation — FEniCSx tutorial

e1_facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, lambda x: electrode1(x, ns))
e2_facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, lambda x: electrode2(x, ns))
wall_facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, lambda x: wall(x, ns))
marked_facets = np.hstack([e1_facets, e1_facets, wall_facets])
marked_values = np.hstack([
    np.full_like(e1_facets, 1),
    np.full_like(e2_facets, 2),
    np.full_like(wall_facets, 3),
])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
ds = ufl.Measure('ds', domain=mesh, subdomain_data=facet_tag)

wall_dofs = dolfinx.fem.locate_dofs_topological(basis, fdim, wall_facets)
bcs = [
    dolfinx.fem.dirichletbc(dolfinx.default_scalar_type(0), wall_dofs, basis),
]

and then writing the kernels as

zsh = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1))
sigma = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(np.eye(3)))
a = ufl.inner(sigma * ufl.grad(u), ufl.grad(v)) * ufl.dx
b = sum([1/zsh * ufl.inner(u, v) * ufl.ds(i) for i in (1,2)])
c = [1/zsh * v * ds(i) for i in (1, 2)]
d_diag = [1/zsh * ds(i) for i in (1, 2)]

I think this implements the A, B, C, * D in the original post. But how to assemble them into one LinearProblem?

Also, if I am doing something that is not going to scale to large meshes well, please adivse.