# Imposing a discontinuity at interface using DG method

**URL:** <https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780>\
**Category:** dolfinx\
**Created:** [July 19, 2023, 7:40pm UTC](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780 "2023-07-19T19:40:17Z")\
**Posts on this page:** 1\
**Showing post:** 2

<div class="post-metadata">

**Author:** ![violetus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/violetus/32/5122_2.png) [@violetus](https://fenicsproject.discourse.group/u/violetus)\
**Post date:** [July 21, 2023, 1:50am UTC](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780/2 "2023-07-21T01:50:51Z")

</div>

Hi, making the following changes made the code work for me -

1. Using `from dolfinx import io` instead of `from dolfinx.io import XDMFFile`
2. Introducing a `tol` at interface like this -

```auto
tol = 1e-8
def marker_interface(x):
    return np.isclose(x[0], 0.5, tol)

```

1. Replacing `dot` with `inner` like -

```auto
# diffusion
F += inner(grad(v), grad(u))*dx - inner(v*n, grad(u))*ds \
   - inner(avg(grad(v)), jump(u, n))*dS - inner(jump(v, n), avg(grad(u)))*dS \
   + gamma/avg(h)*inner(jump(v, n), jump(u, n))*dS

```

and

```auto
F += - inner(grad(v), u*n)*ds + alpha/h*v*u*ds\
   + uD*inner(grad(v), n)*ds - alpha/h*uD*v*ds

```

1. Writing `u` to `XDMFFile` as -

```auto
with io.XDMFFile(MPI.COMM_WORLD, "uh.xdmf", "w") as file:
        file.write_mesh(u.function_space.mesh)
        file.write_function(u, 0)
        file.close()

```

Hence the whole code is -

```auto
from dolfinx import mesh, fem, nls, plot
from mpi4py import MPI
import ufl
from dolfinx import io
from petsc4py import PETSc
from ufl import dx, grad, dot, jump, avg
import numpy as np

msh = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)

V = fem.FunctionSpace(msh, ("DG", 1))

# create mesh tags
tol = 1e-8
def marker_interface(x):
    return np.isclose(x[0], 0.5, tol)

tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
facet_imap = msh.topology.index_map(tdim - 1)
boundary_facets = mesh.exterior_facet_indices(msh.topology)
interface_facets = mesh.locate_entities(msh, tdim - 1, marker_interface)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
indices = np.arange(0, num_facets)
values = np.zeros(indices.shape, dtype=np.intc) # all facets are tagged with zero

values[boundary_facets] = 1
values[interface_facets] = 2

mesh_tags_facets = mesh.meshtags(msh, tdim - 1, indices, values) 

ds = ufl.Measure("ds", domain=msh, subdomain_data=mesh_tags_facets)
dS = ufl.Measure("dS", domain=msh, subdomain_data=mesh_tags_facets)

u = fem.Function(V)
u_n = fem.Function(V)
v = ufl.TestFunction(V)

h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)

# Define parameters
alpha = 10
gamma = 10

# Simulation constants
f = fem.Constant(msh, PETSc.ScalarType(2.0))

# Define variational problem
F = 0

# diffusion
F += inner(grad(v), grad(u))*dx - inner(v*n, grad(u))*ds \
   - inner(avg(grad(v)), jump(u, n))*dS - inner(jump(v, n), avg(grad(u)))*dS \
   + gamma/avg(h)*inner(jump(v, n), jump(u, n))*dS

# source
F += -v*f*dx 

# Dirichlet BC
uD = fem.Function(V)
uD.interpolate(lambda x: np.full(x[0].shape, 0.0))

F += - inner(grad(v), u*n)*ds + alpha/h*v*u*ds\
   + uD*inner(grad(v), n)*ds - alpha/h*uD*v*ds

problem = fem.petsc.NonlinearProblem(F, u)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(u)

with io.XDMFFile(MPI.COMM_WORLD, "uh.xdmf", "w") as file:
        file.write_mesh(u.function_space.mesh)
        file.write_function(u, 0)
        file.close()

```

Visualization in paraview -

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/ae77a0c3a97aa49be60d54e932fa27c54196166e.png)

---

_[View the full topic](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780)._