Imposing a discontinuity at interface using DG method

dolfinx
2

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 -
tol = 1e-8
def marker_interface(x):
    return np.isclose(x[0], 0.5, tol)
  1. Replacing dot with inner like -
# 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

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

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
image1920×862 186 KB