I am solving a variation of the poisson problem with two domains as shown at Defining subdomains for different materials — FEniCSx tutorial. How do I add a discontinuity at the interface of type:
\nabla u_1 \cdot \hat{n} = -\frac{u_0 - u_1}{\kappa_1}
I have tried mixed elements to have u_0 and u_1 separate but my dimensions are all mismatched.
For anyone interested, see solution at Imposing a discontinuity at interface using DG method - #19 by smesc
Hi leshinka,
I’m dealing with the same discontinuity as you are. From the solution proposed on the thread you linked, I cannot directly figure out how to solve our problem. Could you share how you solved it? Specifically, I’m curious what your variational form looks like.
Thanks in advance!
My problem is to prescribe a specific flux condition in a set of internal facets (the internal discontinuity boundary). The way that I chose to solve the entire problem is using Discontinuous Galerkin (DG) with interior penalty. I prescribe Dirichlet and Neumann boundary conditions using Nitsche’s method. Dokkens tutorial shows how to solve a poisson problem with DG/IP and prescribing Dirichlet boundary conditions.
While I have the usual IP for the rest of the internal facets, I set a specific form of flux for the Newton cooling type boundary.
[\![ u ]\!] \overset{\Delta}{=} u(+) \cdot \hat{\mathbf{n}}(+) + u(-) \cdot \hat{\mathbf{n}}(-) = \frac{1}{\kappa_1}\nabla u (-)
Therefore, wherever in my DG formulation I have the jump of test function u, I replace that with \frac{1}{\kappa_1}\nabla u (-). You will also need to edit the coercivity term for the internal discontinuity boundary.
My problem has several Dirichlet and Neumann boundary conditions, so the variational form is long and I don’t have time to simplify it to a MWE. However, you’re welcome to share your MWE and I can help debug.
Thank you for your response! It helped me get to a solution that seems reasonable. However, it seems leaving some terms out, changing some dS to dS(0), and similar, minor changes, also lead to similar, reasonable solutions. Therefore I have difficulty verifying that the current formulation is correct. Could you have a look at it, to see if I’m missing something?
I have a similar boundary condition to yours, where
\theta [[u]] = - k_0 \nabla u \cdot \mathbf{n}(-)
and therefore, the corresponding terms in the variational form are slightly different. The MWE can be found below.
P.S. I’m using FEniCS rather than FEniCSx, due to some dependencies in other packages I use.
from fenics import *
if __name__ == '__main__':
# define global tolerance (i.e. we consider y-tol <= x <= y+tol iff x == y)
tol = 1e-14
# create mesh, function space, trial and test functions
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, 'DG', 1)
u = TrialFunction(V)
v = TestFunction(V)
# subdomains
class Left(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5 + tol
class Right(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= 0.5 + tol
domain_markers = MeshFunction('size_t', mesh, mesh.topology().dim())
subdomain0 = Left()
subdomain1 = Right()
subdomain0.mark(domain_markers, 0)
subdomain1.mark(domain_markers, 1)
class BoundaryLeftTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol) and (0.5 - tol <= x[1] <= 1 + tol)
class BoundaryLeftBottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol) and (0 - tol <= x[1] <= 0.5 + tol)
class BoundaryRight(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, tol)
class BoundaryTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1, tol)
class BoundaryBottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, tol)
class BoundaryMiddle(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.5, tol)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
# heater boundary on left
boundary_heater = BoundaryLeftTop()
boundary_heater.mark(boundary_markers, 0)
# robin boundary on the remaining part
boundary_robin_1 = BoundaryTop()
boundary_robin_1.mark(boundary_markers, 1)
boundary_robin_2 = BoundaryRight()
boundary_robin_2.mark(boundary_markers, 1)
boundary_robin_3 = BoundaryBottom()
boundary_robin_3.mark(boundary_markers, 1)
boundary_robin_4 = BoundaryLeftBottom()
boundary_robin_4.mark(boundary_markers, 1)
boundary_5 = BoundaryMiddle()
boundary_5.mark(boundary_markers, 2)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=domain_markers)
dS = Measure('dS', domain=mesh, subdomain_data=boundary_markers)
# parameters
class K(UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
self.materials = materials
self.k_0, self.k_1 = k_0, k_1
super().__init__(**kwargs)
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
k_0 = 0.5
k_1 = 1
k = K(domain_markers, k_0, k_1, degree=0)
r = Constant(0.1)
s = Constant(1)
g = Constant(1)
# Define normal vector and mesh size
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-')) / 2
gamma = 10.0 # coercivity parameter
theta = 0.5 # middle boundary conductivity
# variational problem
F = k * dot(grad(u), grad(v))*dx + r * (u - s) * v * ds(1) - g * v * ds(0)
# conductivity robin boundary neumann boundary
F += (gamma / h_avg) * dot(jump(u, n), jump(v, n)) * dS # coercivity
F -= dot(avg(grad(u)), jump(v, n)) * dS
F -= dot(jump(u, n), avg(grad(v))) * dS # symmetry
# interface
F -= dot(avg(grad(v)), n('-')) * (dot(grad(u)('-'), n('-')) * (k_0 / theta)) * dS(2)
F += gamma / avg(h) * dot(jump(v, n), n('-')) * (dot(grad(u)('-'), n('-')) * (k_0 / theta)) * dS(2)
a, L = lhs(F), rhs(F)
# compute solution
u = Function(V)
solve(a == L, u)
# save solution to file
vtkfile = File('out/solution.pvd')
vtkfile << u
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