Calculating quantity containing the unit normal on the boundary with Discontinuous Galerkin

Errors
1

Hi all, I am solving the wave equation in a 2d square domain. Inside I have an circle as an obstacle. I want to calculate

stress = dot(grad(u), n)

on the circle boundary inside. When I try to project it I get an error: Here is the code

import platform, sys
python_version=platform.python_version()
from distutils.version import LooseVersion, StrictVersion
from matplotlib import cm, colors
import mshr
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt;
from IPython.display import clear_output, display; import time; import dolfin.common.plotting as fenicsplot 
import time
import os, sys, shutil
from os import path
import numpy as np
import matplotlib as mpl

# solver for the Ax=b system. If not specified, a Newton solver is used
# but GMRES is particularly for non-linear systems
solver = KrylovSolver("gmres")
solver.parameters["relative_tolerance"] = 5e-6
solver.parameters["maximum_iterations"] = 1000
solver.parameters["monitor_convergence"] = False


a_c = 1 #  radius
a_i = 0.5
x_c = 0.0
y_c = 0.0
x_s = 1.0 # point source coord.
y_s = 1.0

X1, Y1 = -2.0,-2.0;
X3, Y3 =  2.0, 2.0;

# Geometry of the background/soil domain
domain = mshr.Rectangle(dolfin.Point(X1, Y1), dolfin.Point(X3, Y3))

# Geoemtry for the obstacle 
domain_obstacle = mshr.Circle(dolfin.Point(x_c, y_c), a_c)
# Set subdomain for the obstacle
domain.set_subdomain(1, domain_obstacle)

n_elem = 2**6
mesh = mshr.generate_mesh(domain, n_elem)

# Time variables
T = 100
Ndt = 2000
dt = T / Ndt
#c = 1
tol = 1E-10

# Define classes for the different boundaries in the domain
class Obstacle_out(SubDomain): # also used for inclusion
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def inside(self, x, on_boundary):
        tol = 1E-2   # tolerance for coordinate comparisons: here for less it does not find the boundary
        return near(x[0]*x[0] + x[1]*x[1], a_c*a_c, tol)

class Obstacle_in(SubDomain):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def inside(self, x, on_boundary):
        tol = 1E-2   # tolerance for coordinate comparisons: here for less it does not find the boundary
        return near(x[0]*x[0] + x[1]*x[1], a_i*a_i, tol)        

class Left(SubDomain):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def inside(self, x, on_boundary):
        return near(x[0], -2.0, tol)

class Right(SubDomain):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def inside(self, x, on_boundary):
        return near(x[0], 2.0, tol)

class Top(SubDomain):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def inside(self, x, on_boundary):
        return near(x[1], 2.0, tol) and on_boundary

class Bottom(SubDomain):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def inside(self, x, on_boundary):
        return near(x[1], -2.0, tol) and on_boundary

markers = MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
facetmarkers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1,0)

# Initialize and mark subdomains
facetmarkers.set_all(0)

top    = Top(markers)
bottom = Bottom(markers)
left   = Left(markers)
right  = Right(markers)
obstacle_out = Obstacle_out(markers)
obstacle_in = Obstacle_in(markers)

# assigned markers: 0 for medium, 1 for obstacle

top.mark(facetmarkers,0)
bottom.mark(facetmarkers,1) #Dirichlet boundary
left.mark(facetmarkers,2)
right.mark(facetmarkers,3)
obstacle_out.mark(facetmarkers,4)
obstacle_in.mark(facetmarkers,5) #Dirichlet boundary

dx = Measure('dx', domain=mesh, subdomain_data=markers) 
ds = Measure('ds', domain=mesh, subdomain_data=facetmarkers) 
dS = Measure('dS', domain=mesh, subdomain_data=facetmarkers)

rho = Constant(1.0)
mu = Constant(1.0)
c = sqrt(mu/rho)

# Define the function space
V = FunctionSpace(mesh, "DG", 1) 
Vc = FunctionSpace(mesh, "CG", 1)

# Dirichlet at the bottom
uD = Expression('0.0', degree=2)

# initial condition, previous and current solution
g_ini = Constant(tol)
u1= interpolate(g_ini, V)
u0= interpolate(g_ini, V)

# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)

n = FacetNormal(mesh)
h = CellDiameter(mesh)

alpha = 4
gamma = 8
K1 = 4
K2 = 2

# domain integrals
a1 = inner(u,v)*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx \
     - dt*dt*c*c*dot(avg(grad(v)), jump(u, n))*dS(0) \
     - dt*dt*c*c*dot(jump(v, n), avg(grad(u)))*dS(0) \
     + dt*dt*c*c*alpha/avg(h)*dot(jump(v, n), jump(u, n))*dS(0)

# Dirichlet boundary Nitsche terms
a2 = - dt*dt*c*c*dot(grad(v), u*n)*ds(1) \
     - dt*dt*c*c*dot(v*n, grad(u))*ds(1) \
     + dt*dt*c*c*(gamma/h)*v*u*ds(1)

# Interface terms     
a3 = - dt*dt*c*c*dot(avg(grad(v)), n('-'))*(u('-')*(K1/K2-1))*dS(4)\
     + dt*dt*c*c*alpha/avg(h)*dot(jump(v,n),n('-'))*(u('-')*(K1/K2-1))*dS(4)

# symmetry
a4 = - dt*dt*c*c*dot(avg(grad(v)), jump(u, n))*dS(4)
# coercivity
a5 = + dt*dt*c*c*alpha/avg(h)*dot(jump(v, n), jump(u, n))*dS(4)

a  = a1+a2+a3+a4+a5

L  = + 2*inner(u1,v)*dx-inner(u0,v)*dx  \
     - dt*dt*c*c*uD*dot(grad(v), n)*ds(1) + dt*dt*c*c*(gamma/h)*uD*v*ds(1) \

u = Function(V)

stepcounter = 0;
t = 0;
while t <= T:    
    A, b = assemble_system(a, L)
    deltasource = PointSource(V, Point(x_s, y_s), t*exp(-t))
    deltasource.apply(b)
    solver.solve(A, u.vector(), b)
    u0.assign(u1)
    u1.assign(u)
   
    stress = dot(grad(u), n) 
    #stress_project = project(stress, Vc) # this gives me error
    plot(u)
    plt.pause(0.002)
    stepcounter += 1;    
    t += dt
2

This doesn’t make sense, as a standard projection is in L^2(\Omega), while the normal vector only lives on facets.
Also note that Vc should be discontinuous, as the gradient of a continuous space is discontinuous.

See for instance: Obtain velocity normal to boundary - #2 by dokken