Weak imposition of boundary condition on interior nodes

variational formulation
1

Hello,

This is a follow-up question to: Using different time steps in different subdomains

I have different subdomains and need to apply a Dirichlet boundary condition on all the interior nodes of one of these subdomains.

I figured out how to do this using strong imposition by checking for all DOFs in the respective subdomain. Is there a way to do this using weak imposition and DG-elements?

I can implement boundary conditions using the integration measure ds on external boundaries (as it is shown in the tutorial). Trying the same with dS for the interior case did not work. I tried for example dS(2) for the subdomain that is marked as 2.

Thank you in advance!

2

To impose boundary conditions weakly with DG finite elements, you would have to think carefully about the variational formulation, as a DG-formulation it-self uses weak imposition of continuity (a Nitsche-like method).

Without an minimal code example illustrating your problem, it is hard to give you any guidance.

3

I try to simplify my problem in the following example. Here I solve the heat equation with two different materials with a discontinuous method using the following domain:

image

Here is the code:

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx import fem, mesh, io, plot
from ufl import dot, grad, jump, avg, FacetNormal, Circumradius, dx, ds, dS

# Define temporal parameters
t = 0 # Start time
T = 2.0 # Final time
num_steps = 61     
dt = T / num_steps # time step size

# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])], [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("DG", 1))

Q = fem.FunctionSpace(domain, ("DG", 0))

def Omega_1(x):
    return x[1] >= 0.0

def Omega_2(x):
    return x[1] <= 0.0

kappa = fem.Function(Q)
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)
cells_2 = mesh.locate_entities(domain, domain.topology.dim, Omega_2)

kappa.x.array[cells_1] = np.full_like(cells_1, 1, dtype=ScalarType)
kappa.x.array[cells_2] = np.full_like(cells_2, 0.01, dtype=ScalarType)

# Create initial condition
def initial_condition(x, a=5):
    return np.exp(-a*(x[0]**2+x[1]**2))
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

xdmf = io.XDMFFile(domain.comm, "output.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)

import ufl
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
uD = fem.Constant(domain, PETSc.ScalarType(0))
n = FacetNormal(domain)

# bilinear form
a = u * v * dx + dt * dot(grad(u), grad(kappa * v)) * dx
a += dt * (dot(jump(u, n), avg(grad(kappa * v))) - dot(avg(grad(kappa * u)), jump(v, n))) * dS
a += dt * (dot(u * n, grad(kappa * v)) - dot(grad(kappa * u), v * n)) * ds
a += dt * dot(jump(u, n), jump(v, n)) * dS
a += dt * dot(u * n, v * n) * ds

# liniear form
L = (u_n + dt * f) * v * dx
L += dt * dot(uD * n, grad(kappa * v)) * ds
L += dt * uD * v * ds

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = fem.petsc.assemble_matrix(bilinear_form)
A.assemble()
b = fem.petsc.create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    fem.petsc.assemble_vector(b, linear_form)

    # Solve linear problem
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array

    # Write solution to file
    xdmf.write_function(uh, t)

xdmf.close()

I use a weak imposition of the Dirichlet boundary which is zero on the external boundary of the domain. What I need to know is how to modify / add a term to my variational formulation to impose a boundary condition on interior nodes, in this case subdomain 2 to fixate the solution there.

For the variational formulation I use an interior penalty method with all parameters set to 1 for simplicity.