I am pretty new to FEniCS and to FEA, so I skipped fenics
and went directly to dolfinx
. However, perhaps I can point out a couple of things. Perhaps you can split the variational form from the boundary conditions into separate lines (for readability and ease of debug).
From the tutorial at https://github.com/jorgensd/dolfinx-tutorial/blob/dokken/dg_tutorial/chapter1/dg.ipynb, I think your variational formulation should be
k * dot(grad(u), grad(v))*dx - dot(n, grad(u)) * v * ds
Otherwise, I am not entire sure that you can set up the boundary conditions in DG the way you’ve done. The tutorial at https://github.com/jorgensd/dolfinx-tutorial/blob/dokken/dg_tutorial/chapter1/nitsche.ipynb shows how to set up boundary conditions using Nitsche’s method.
If you’re using interior penalty then you’re gonna have to mark the rest of the interior facets (and not just the middle boundary). I don’t know how to do this in fenics
, but here is how I did it in dolfinx
from already labelled exterior boundary and middle boundary (really shared boundary between two marked subdomains):
comm = MPI.COMM_WORLD
with io.XDMFFile(comm, tria_meshfile, "r") as infile3:
domain = infile3.read_mesh(cpp.mesh.GhostMode.none, 'Grid')
ct = infile3.read_meshtags(domain, name="Grid")
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(tdim, fdim)
ft_imap = domain.topology.index_map(fdim)
num_facets = ft_imap.size_local + ft_imap.num_ghosts
indices = np.arange(0, num_facets)
values = np.zeros(indices.shape, dtype=np.intc) # all facets are tagged with zero
with io.XDMFFile(comm, line_meshfile, "r") as infile2:
ft = infile2.read_meshtags(domain, name="Grid")
values[ft.indices] = ft.values
meshtags = mesh.meshtags(domain, fdim, indices, values)
domaintags = mesh.meshtags(domain, domain.topology.dim, ct.indices, ct.values)
If you have anaconda, I think you can install dolfinx quite easily in a custom environment according to instructions at GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment. Then you can customize the solution given by smesc at Imposing a discontinuity at interface using DG method - #19 by smesc