I am in the process of translating a working code from FEniCS to FEniCSx, but I am stuck at the boundary conditions part.
First, the mesh comes from a file.geo created with GMSH GUI:
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1e-2, 1e-2, 1};
Physical Volume("the_wire", 13) = {1};
Physical Surface("left_end", 14) = {5};
Physical Surface("right_end", 15) = {6};
It’s just a box that represents a wire. The idea is to solve Ohm’s law: we apply a known current to the wire, and “measure” its voltage, this should yield the resistance of the wire. More precisely, the current is injected via the left end of the wire (surface 5), and leaves the wire at the right end (surface 6). The voltage is read between these two surfaces, too.
I now show you my code in 2 parts, the first part works, I believe, and is just there as to make a MWE.
1st part:
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("file.msh", MPI.COMM_WORLD, gdim=3)
# Now I want to convert the mesh file to XDMF using meshio.
proc = MPI.COMM_WORLD.rank
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
if proc == 0:
# Read in mesh
msh = meshio.read("file.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
tetra_mesh = create_mesh(msh, "tetra", prune_z=False)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mt.xdmf", triangle_mesh)
# Now read the mesh files produced above.
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mesh, name="Grid")
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner, as_tensor, inv)
import numpy as np
# Define function space
V = FunctionSpace(mesh, ("CG", 1))
volt = Function(V)
u = TestFunction(V)
rho_tensor = as_tensor([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
sigma_tensor = inv(rho_tensor)
# Reminder of the mesh:
# Physical Volume("the_wire", 13) = {1};
# Physical Surface("left_end", 14) = {5};
# Physical Surface("right_end", 15) = {6};
inj_current_surface = 5
vanish_voltage_surface = 6 # This corresponds to the surface the current leaves the material.
V_current_contact_out = 0.0 # Voltage value of the surface where the current leaves the material.
reading_voltage_surface_0 = 5
reading_voltage_surface_1 = 6
Now the 2nd part, where I am stuck.
# We define the boundary conditions
u_bc = Function(V)
left_facets = ft.find(inj_current_surface)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bcs = [dirichletbc(ScalarType(1), left_dofs, V)]
Here I have the first doubt, I want to apply the Dirichlet boundary condition of voltage = 0 on the left end of the wire. In FEniCS, I have the code:
bc0 = DirichletBC(V, V_current_contact_out, boundary_parts, vanish_voltage_surface)
bcs = [bc0]
J = 2e-3 / assemble(1 * ds(inj_current_surface, domain=mesh))
F = dot(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(inj_current_surface)
solve(F == 0, volt, bcs)
That’s the part I am stuck at translating. Any help is greatly appreciated.