Is there a tutorial on electric field for a parallel plate capacitor using Fenicsx?

General
1

I want to learn to model electric field around electrostatic charged bodies using Fenicsx and looking for tutorials on parallel capacitors and similar problems.

2 
import numpy as np
import ufl
from dolfinx.fem import functionspace, Function, dirichletbc
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from mpi4py import MPI
import matplotlib.pyplot as plt
from dolfinx import default_scalar_type, fem, mesh, io
from petsc4py import PETSc

# Create mesh and define function space
nx, ny = 50, 50

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])],
                               [nx, ny], mesh.CellType.triangle)
V = fem.functionspace(domain, ("CG", 1))

f = fem.Constant(domain, default_scalar_type(0.0))

# Define boundary conditions
def boundary_top(x):
    return np.isclose(x[1], 2)

def boundary_bottom(x):
    return np.isclose(x[1], -2)

# Set Dirichlet boundary conditions
fdim = domain.topology.dim - 1
facet_top = mesh.locate_entities_boundary(domain, fdim, boundary_top)


# Define a function to set the inflow profile
def sin_function(x):
    return np.sin(np.pi * x[1])

u_top = fem.Function(V)
u_top.interpolate(sin_function)
dofs_top = fem.locate_dofs_topological(V, fdim, facet_top)

bc_top = fem.dirichletbc(u_top, dofs_top, V)

facet_bottom = mesh.locate_entities_boundary(domain, fdim, boundary_bottom)
bc_bottom = fem.dirichletbc(PETSc.ScalarType(-1), fem.locate_dofs_topological(V, fdim, facet_bottom), V)


# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
 
 
# Compute solution
A_z = fem.Function(V)
problem = LinearProblem(a, L, u=A_z, bcs=[bc_top, bc_bottom], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# Save solution to file
with io.XDMFFile(domain.comm, "capacitor_potential.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(uh)

# Compute electric field
gdim = domain.geometry.dim
E_field = -ufl.grad(uh)
V_E = fem.functionspace(domain, ("CG", 2, (gdim, )))
E = fem.Function(V_E)
E_expr = fem.Expression(E_field, V_E.element.interpolation_points())
E.interpolate(E_expr)


# Export solution
with io.VTXWriter(domain.comm, "capacitor_field.bp", E, engine="BP4") as vtx:
    vtx.write(0.0)

I have written a small code, at the top I want a sinusoidal distribution of potential. But I am getting some error. Also, how can I expand the domain to include space around the plate?

3

Since no tutorial was available for solving basic electrostatics problem using Fenics. I tried to try it for solving capacitor.

import numpy as np
import ufl
import gmsh
import dolfinx
from petsc4py import PETSc
from dolfinx import fem, mesh, io
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from dolfinx.io import gmshio
from ufl import div, grad, inner, dx, Measure

# Initialize parameters
L, H = 5, 5  # Domain size
plate_length = 1
plate_height = 0.1
V0 = 100.0  # Voltage on upper plate
model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 2

# Initialize gmsh
gmsh.initialize()

if mesh_comm.rank == model_rank:
    # Create background domain
    background = gmsh.model.occ.addRectangle(0, 0, 0, L, H)
    
    # Create plates
    upper_plate = gmsh.model.occ.addRectangle(2, 2, 0, plate_length, plate_height)
    lower_plate = gmsh.model.occ.addRectangle(2, 3, 0, plate_length, plate_height)
    gmsh.model.occ.synchronize()
    
    # Fragment all geometries together
    all_surfaces = [(2, upper_plate), (2, lower_plate)]
    whole_domain = gmsh.model.occ.fragment([(2, background)], all_surfaces)
    gmsh.model.occ.synchronize()
    
    # Create physical markers
    upper_plate_marker, lower_plate_marker = 1, 2
    background_surfaces = []
    other_surfaces = []
    
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        
        # Identify upper plate
        if 2 <= com[0] <= 3 and 2 <= com[1] <= 2.1:
            gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=1)
            gmsh.model.setPhysicalName(domain[0], 1, "upper_plate")
            other_surfaces.append(domain)
            
        # Identify lower plate
        elif 2 <= com[0] <= 3 and 3 <= com[1] <= 3.1:
            gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=2)
            gmsh.model.setPhysicalName(domain[0], 2, "lower_plate")
            other_surfaces.append(domain)
            
        # Background is everything else
        else:
            background_surfaces.append(domain[1])
    
    # Add marker for background
    gmsh.model.addPhysicalGroup(2, background_surfaces, tag=0)
    gmsh.model.setPhysicalName(2, 0, "background")
    
    # Create refined mesh around plates
    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(other_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
    
    # Add threshold field for mesh refinement
    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", plate_height/3)
    gmsh.model.mesh.field.setNumber(2, "LcMax", plate_height*6)
    gmsh.model.mesh.field.setNumber(2, "DistMin", plate_height*4)
    gmsh.model.mesh.field.setNumber(2, "DistMax", plate_height*10)
    
    gmsh.model.mesh.field.add("Min", 5)
    gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
    gmsh.model.mesh.field.setAsBackgroundMesh(5)
    
    # Generate mesh
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")
    gmsh.write("capacitor.msh")

# Read mesh
mesh_obj, cell_tags,facet_tags  = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()

with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_obj)
    xdmf.write_meshtags(cell_tags, mesh_obj.geometry)

I was able to create mesh. But not able to apply dirichletbc

# Create function space (using 2nd order elements for better accuracy)
V = fem.functionspace(mesh_obj, ("CG", 1))

# Define boundary conditions
upper_plate_dofs = fem.locate_dofs_topological(V, gdim, cell_tags.find(1))
lower_plate_dofs = fem.locate_dofs_topological(V, gdim, cell_tage.find(2))

# Set Dirichlet boundary conditions
bc_upper = fem.dirichletbc(PETSc.ScalarType(V0), upper_plate_dofs, V)
bc_lower = fem.dirichletbc(PETSc.ScalarType(-V0), lower_plate_dofs, V)
4

The mesh generation is quite involved, so I didn’t look through all details.
Still, your definition of boundary conditions is certainly wrong. You need to tag entities of dimension gdim - 1 to apply boundary conditions on them.