Hi all,
This is a simpler version of the problem stated in Imposing a discontinuity in the solution, on an internal interface
I’d like (in dolfinx) to impose a discontinuity at the interface between two domains.
I understood that MeshView would be the best option since we would have CG elements everywhere except at the interface. However it’s not yet implemented in FEniCSx
I’d like to give the DG method a go!
This is a simple poisson problem and I’d like to enforce the following equation at the interface x=0.5:
\frac{u(-)}{K(-)} = \frac{u(+)}{K(+)}
\nabla u(-) - \nabla u(+) = 0
Bonus: I would be curious to also have the solution for \frac{u(-)}{K(-)} = \frac{u(+)^2}{K(+)^2}
Where u is my solution and K are constants.
I’ve got the DG method working for the simple poisson problem, I’m just struggling to find how to enforce the interface condition.
This is my MWE:
from dolfinx import mesh, fem, nls, plot
from mpi4py import MPI
import ufl
from dolfinx.io import XDMFFile
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
def marker_interface(x):
return np.isclose(x[0], 0.5)
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 += dot(grad(v), grad(u))*dx - dot(v*n, grad(u))*ds \
- dot(avg(grad(v)), jump(u, n))*dS - dot(jump(v, n), avg(grad(u)))*dS \
+ gamma/avg(h)*dot(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 += - dot(grad(v), u*n)*ds + alpha/h*v*u*ds\
+ uD*dot(grad(v), n)*ds - alpha/h*uD*v*ds
# jump at interface
# F += ...... * dS(2)
problem = fem.petsc.NonlinearProblem(F, u)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(u)
xdmf_file = XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w")
xdmf_file.write_mesh(msh)
xdmf_file.write_function(u)
EDIT: With the above MWE, I obtain the following solution (without the discontinuity at x=0.5)
Thanks in advance for your help!
