# Poisson equation with two subdomains and Newton cooling type discontinuity at the interface

**URL:** <https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673>\
**Category:** General\
**Created:** [November 3, 2022, 9:59pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673 "2022-11-03T21:59:15Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [November 3, 2022, 9:59pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673/1 "2022-11-03T21:59:15Z")

</div>

I am solving a variation of the poisson problem with two domains as shown at [Defining subdomains for different materials — FEniCSx tutorial](https://jorgensd.github.io/dolfinx-tutorial/chapter3/subdomains.html). 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.

---

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [October 31, 2023, 11:31pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673/2 "2023-10-31T23:31:13Z")

</div>

For anyone interested, see solution at [Imposing a discontinuity at interface using DG method - #19 by smesc](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780/19)

---

<div class="post-metadata">

**Author:** ![basgieling](https://avatars.discourse-cdn.com/v4/letter/b/e36b37/32.png) [@basgieling](https://fenicsproject.discourse.group/u/basgieling)\
**Post date:** [November 10, 2023, 10:34am UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673/3 "2023-11-10T10:34:09Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [November 10, 2023, 2:53pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673/4 "2023-11-10T14:53:59Z")

</div>

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](https://github.com/jorgensd/dolfinx-tutorial/blob/dokken/dg_tutorial/chapter1/dg.ipynb) 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.

---

<div class="post-metadata">

**Author:** ![basgieling](https://avatars.discourse-cdn.com/v4/letter/b/e36b37/32.png) [@basgieling](https://fenicsproject.discourse.group/u/basgieling)\
**Post date:** [November 13, 2023, 1:16pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673/5 "2023-11-13T13:16:12Z")

</div>

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.

```auto
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

```

---

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [November 13, 2023, 2:22pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-two-subdomains-and-newton-cooling-type-discontinuity-at-the-interface/9673/6 "2023-11-13T14:22:22Z")

</div>

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](https://github.com/jorgensd/dolfinx-tutorial/blob/dokken/dg_tutorial/chapter1/dg.ipynb), I think your variational formulation should be

```auto
 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](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):

```auto
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](https://github.com/FEniCS/dolfinx). Then you can customize the solution given by smesc at [Imposing a discontinuity at interface using DG method - #19 by smesc](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780/19)