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

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.

Thanks in advance!

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]

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",
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):
        mask = CellFunctionSizet(mesh, 1)
        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)]):
              mask.set_value(c.index(), 1)
              mask.set_value(c.index(), 0)
        self.mask = mask
    def eval_cell(self, values, x, ufc_cell):
        if self.mask.array()[ufc_cell.index] == 1:
            values[0] = 1
            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)
mask = interpolate(jB, V) 

F1 = interpolate(Constant(1.0), V)     
F2 = project(jump(F1)*mask, V) 
written 6 months ago by Dan Sweeney  

1 Answer

6 months ago by
I also found something like this that may work
Vs = FunctionSpace(mesh, 'DG', 0) 
v0 = TestFunction(Vs)
mask_fxn = project(Constant(1.0), Vs)
mask = Function(Vs, assemble(mask_fxn('+')*v0('+')*dS(7)))​
Please login to add an answer/comment or follow this question.

Similar posts:
Search »