Calculating electrostatic force through virtual work

variational formulation
1

Hi,

I want to compute the electrostatic force through virtual work method. This involves calculating the sensitivity of energy to a virtual displacement of the object the force is calculated on.

I am not knowledgeable in these things but it seems that shape derivatives could help. After reading a lot of forum posts on the subject, I am still unable to piece things together to solve the entire problem. Specifically, I do not feel sure how to define the displacement field for this case which defines the derivative and then of course its implementation.

Below I am sharing an MWE involving two circles/cylinders in electrostatic field. I need the force on only one of the two circles so only that should be subject to virtual displacement. I will sincerely appreciate if someone could complete the code below or refer me to examples that realize something similar.

import dolfinx, ufl, mpi4py, petsc4py, gmsh, typing
from dolfinx.io import gmshio
from dolfinx.fem.petsc import LinearProblem
import numpy as np

y0 = 0 # sensitivity parameter

#%% geometry 
gdim = 2
model_rank = 0 # everything is on one core

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
occ = gmsh.model.occ

c1 = occ.addDisk(0, 1+y0, 0, 0.5, 0.5)
c2 = occ.addDisk(0, -1, 0, 0.5, 0.5)
r1 = occ.addRectangle(-3, -3, 0, 6, 6)

frag, _ = occ.fragment([(gdim, c1)], [(gdim, r1), (gdim, c2)])
occ.synchronize()

# number all domains
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)
# gmsh.fltk.run()
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
gmsh.finalize()

dom_indices = (1, 2, 3) # top circle, bottom circle, surrounding domain

# assign permittivity to circles
Omega_dg = dolfinx.fem.functionspace(mesh, ("DG", 0))
epsr = dolfinx.fem.Function(Omega_dg)
epsr.x.array[:] = 1
for tag in dom_indices[:2]:
    cells = ct.find(tag)
    epsr.x.array[cells] = 9

Omega = dolfinx.fem.functionspace(mesh, ("CG", 2))

u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)

e = epsr*ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx - ufl.dot(dolfinx.fem.Constant(mesh, 0.0), v)*ufl.dx

a = ufl.lhs(e)
L = ufl.rhs(e)

# boundary conditions
on_left = lambda x: np.isclose(x[0], -3)
dofs_l = dolfinx.fem.locate_dofs_geometrical(Omega, on_left)
dbc_l = dolfinx.fem.dirichletbc(1.0, dofs_l, Omega)
on_right = lambda x: np.isclose(x[0], 3)
dofs_r = dolfinx.fem.locate_dofs_geometrical(Omega, on_right)
dbc_r = dolfinx.fem.dirichletbc(0., dofs_r, Omega)

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem = LinearProblem(a, L, bcs=[dbc_l, dbc_r], petsc_options=petsc_options)

uh = problem.solve()
e_sol = -ufl.grad(uh) # electric field

J = 1/2*epsr*ufl.dot(e_sol, e_sol)*ufl.dx # energy density whose sensitivity is to be calculated

# How to calculate dJ/dy0?
2

I continue to toy with this problem and was trying to work out the code of listing 3 from this paper. I am getting an error at the very end (line 11 of listing 3) and I do not figure out why. Perhaps the syntax has changed? I will appreciate if someone knowledgeable could help me out with the problem.

Below is the MWE:

import dolfinx, ufl, mpi4py
from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np

mesh = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 10, 10)

Omega = dolfinx.fem.functionspace(mesh, ("CG",2))

u = dolfinx.fem.Function(Omega)
v = ufl.TestFunction(Omega)
uh = dolfinx.fem.Function(Omega)

R = ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx + ufl.dot(dolfinx.fem.Constant(mesh, 0.0), v)*ufl.dx

u_dbc1 = dolfinx.fem.Function(Omega)
u_dbc1.x.array[:] = 1
u_dbc2 = dolfinx.fem.Function(Omega)
u_dbc2.x.array[:] = 0

left_boundary = lambda x: np.isclose(x[0], 0)
right_boundary = lambda x: np.isclose(x[0], 1)

left_entities = dolfinx.mesh.locate_entities_boundary(mesh, 1, left_boundary)
left_boundary_dofs = dolfinx.fem.locate_dofs_topological(Omega, 1, left_entities)
right_entities = dolfinx.mesh.locate_entities_boundary(mesh, 1, right_boundary)
right_boundary_dofs = dolfinx.fem.locate_dofs_topological(Omega, 1, right_entities)
dbc = [
    dolfinx.fem.dirichletbc(u_dbc1, left_boundary_dofs),
    dolfinx.fem.dirichletbc(u_dbc2, right_boundary_dofs)]

# fwd_problem = LinearProblem(ufl.lhs(R), ufl.rhs(R), bcs=dbc)
# uh = fwd_problem.solve()
problem_fwd = NonlinearProblem(R, u, dbc)
solver_fwd = NewtonSolver(mpi4py.MPI.COMM_WORLD, problem_fwd)
solver_fwd.rtol = 1e-16
solver_fwd.atol = 1e-10
solver_fwd.solve(u)

J = 1/2*ufl.dot(ufl.grad(u), ufl.grad(u))*ufl.dx
J_form = dolfinx.fem.form(J)
J_val = dolfinx.fem.assemble_scalar(J_form)

# Implementing code snippet from https://doi.org/10.1007/s00158-019-02281-z
dJ_du = ufl.derivative(J, u, v)
dR_du = ufl.derivative(R, u)
lhs = ufl.adjoint(dR_du)
rhs = -dJ_du

problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=[], 
                                          petsc_options={"ksp_type": "preonly", 
                                                         "pc_type": "lu"
                                                         })
lmbda = problem_adj.solve()
x, y = X = ufl.SpatialCoordinate(mesh)
L = ufl.replace(R, {v: lmbda})+J
dJ = ufl.derivative(L, X)

The very last line generates the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[45], line 57
     55 x, y = X = ufl.SpatialCoordinate(mesh)
     56 L = ufl.replace(R, {v: lmbda})+J
---> 57 dJ = ufl.derivative(L, X)

File ~/miniconda3/envs/fenics/lib/python3.12/site-packages/ufl/formoperators.py:349, in derivative(form, coefficient, argument, coefficient_derivatives)
    344     else:
    345         raise NotImplementedError(
    346             "Action derivative not supported when the left argument is not a 1-form."
    347         )
--> 349 coefficients, arguments = _handle_derivative_arguments(form, coefficient, argument)
    350 if coefficient_derivatives is None:
    351     coefficient_derivatives = ExprMapping()

File ~/miniconda3/envs/fenics/lib/python3.12/site-packages/ufl/formoperators.py:213, in _handle_derivative_arguments(form, coefficient, argument)
    208 if argument is None:
    209     # Try to create argument if not provided
    210     if not all(
    211         isinstance(c, (Coefficient, Cofunction, BaseFormOperator)) for c in coefficients
    212     ):
--> 213         raise ValueError(
    214             "Can only create arguments automatically for non-indexed coefficients."
    215         )
    217     # Get existing arguments from form and position the new one
    218     # with the next argument number
    219     if isinstance(form, Form):

ValueError: Can only create arguments automatically for non-indexed coefficients.
3

Consider the following minor adaptation:

Q = dolfinx.fem.functionspace(mesh, mesh.ufl_domain().ufl_coordinate_element())
q = ufl.TestFunction(Q)
dJ = ufl.derivative(L, X, q)
4

Thank you! Could you say a bit to help me understand this step which is very new to me:

Q = dolfinx.fem.functionspace(mesh, mesh.ufl_domain().ufl_coordinate_element())

How to know when we have to define the functionspace this way?

5

You want to differentiate with respect to the nodes of your mesh. These nodes are part of some function space, usually a P-th order continuous Lagrange space.

This used to pull back and push forward values from a physical element to the reference element