I would like to integrate discontinuous functions over a line for a given 2D domain (In practice I am working in 3D so I want to compute the surface integral but convenience I present the 2D case). Here one has to mark the internal edges and then to integrate with the appropriate labels. In the following example I generate a DG function where u = 0 for x[0] < 0.5 and u = 1 for x[0] >=0.5 as one can see from the attached picture
Then integrating u(‘-’) over the line x=0.5 the values should be 1.
In parallel there is some issue when I am using the command
mpirun -np N python example.py
The value of the computed integral is incorrect when N=9 or N =11. Here a minimal example is presented
from dolfin import *
parameters["ghost_mode"] = "shared_facet"
class DiscontinuousOverlLine(UserExpression):
def __init__(self, **kwargs):
try:
super().__init__(**kwargs)
except:
print('super().init failed')
def eval(self, values, x):
values[0] = 1
if x[0] <0.5:
values[0] = 0
return
def value_shape(self):
return ()
M = 2**3
mesh = UnitSquareMesh(M,M)
interior_surface = CompiledSubDomain('near(x[0], 0.5, DOLFIN_EPS)')
interior_surface_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
interior_surface_marker.set_all(0)
interior_surface.mark(interior_surface_marker,1)
V = FunctionSpace(mesh, 'DG', 0)
ut = Function(V)
ut.interpolate(DiscontinuousOverlLine() )
utF = File('ut.pvd')
utF << ut
L =ut('-')*dS(1, subdomain_data=interior_surface_marker)
print('L=', assemble(L))
L =Constant(1)*dS(1, subdomain_data=interior_surface_marker, domain=mesh)
print('L=', assemble(L))