### How do you restrict a DG function to only an internal boundary?

307
views
-1
7 months ago by
I want to do a calculation on an internal interface (labelled 7) using a jump condition between two subdomains using a DG formulation. I then then want to project that to a Function where the values on the interface are the only non-zero values so that I can store it and easily view it in ParaView or something. I'm trying to make sense of previous posts (here and here) toward this end. The code below is what I've come up with, but I'm having trouble using the boundary labels to identify the vertices on the internal interface. I'm wondering if it has to do with DOF mapping in the DG formulation. My mesh files were generated in GMSH and converted with dolfin-convert.

from fenics import *

ofilename = 'simpleModel.xml'
mesh = Mesh(ofilename)
subdomains = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])

V = FunctionSpace(mesh, 'DG', 2)

boundary_verticies = np.where(boundaries.array() == 7)[0]
F1 = interpolate(Constant(1.0), V)
F2 = project(jump(F1), V)
f1 = F1.vector().get_local()
f2 = F2.vector().get_local()
F1.vector()[boundary_vertices] = F2.vector()[boundary_vertices]
plot(F1)​

File attached: simpleModel.xml (238.55 KB)
File attached: simpleModel_facet_region.xml (116.05 KB)

File attached: simpleModel_physical_region.xml (75.57 KB)

Community: FEniCS Project
It's hard to help when the code is not executable due to missing files. See "minimal working example", https://www.allanswered.com/post/emre/read-before-posting-how-do-i-get-by-my-question-answered/
written 6 months ago by Jan Blechta
My apologies--I have added the missing mesh files with a simplified mesh and updated the text. Any advice or suggestions you have would be greatly appreciated!
written 6 months ago by Dan Sweeney
The issue I'm running into in my larger application code is that when I enforce Dirichlet BC's, my jumps become very large around these boundaries. This produces NaNs when I use them as exponents, as I want to do. Ideally, I want to restrict my jump function (or the non-zero values of a projection/interpolation) to the internal domain. Is this possible to do with the mesh from the example I gave above? I've constructed a crude masking function for the boundary in my example code, but this seems very inelegant and I'm wondering if there is a better way to do something similar.

from fenics import *
import numpy as np

class JumpBoundary(Expression):
def __init__(self, mesh, boundaries, boundID, **kwargs):
bf = FacetFunction("uint", mesh, 0)
boundary = np.where(boundaries.array()==boundID)[0]
bf.array()[boundary] = 1
for c in cells(mesh):
if any([c.entities(1)[i] in boundary for i in range(3)]):
else:

def eval_cell(self, values, x, ufc_cell):
values[0] = 1
else:
values[0] = 0

ofilename = 'simpleModel.xml'
mesh = Mesh(ofilename)
subdomains = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])

V = FunctionSpace(mesh, 'DG', 2)

jB = JumpBoundary(mesh, boundaries, 7, degree=0)

F1 = interpolate(Constant(1.0), V)
plot(F2)
written 6 months ago by Dan Sweeney

Vs = FunctionSpace(mesh, 'DG', 0)
mask = Function(Vs, assemble(mask_fxn('+')*v0('+')*dS(7)))​