Point integrals

Question

I'd like to be able to do point integrals like I would boundary or domain, but I can't seem to get it working. In the Demo below, I only get '0.0' for the point integrals (everything else works fine). Any ideas?

Thanks

from dolfin import *

mesh = RectangleMesh(0,0,10,10, 10,10)
V1 = FunctionSpace(mesh, "Lagrange", 1)

u = interpolate(Expression('x[0]+10'),V1)
left = CompiledSubDomain("(std::abs(x[0]) < DOLFIN_EPS) && on_boundary")
corner = CompiledSubDomain("(std::abs(x[0]) < DOLFIN_EPS) && (std::abs(x[1]) < DOLFIN_EPS) && on_boundary")

boundaries = MeshFunction('size_t', mesh, 1) # Gets object of boundaries
boundaries.set_all(0)
left.mark(boundaries,1)
ds = Measure('ds')[boundaries]

points = MeshFunction('size_t', mesh, 0)
points.set_all(0)
corner.mark(points,1)
dp = Measure('point')[points]


print assemble(u*dx)
print assemble(u*ds(1)), assemble(u*ds(0))
print assemble(u*dp(1)), assemble(u*dp(0))

plot(u, interactive=True)

Answer 1

Point integrals does only work within the PointIntegralSolver. It would be nice to get an assemble over either the mesh vertices and/or any given points as you show in your example. In the meanwhile it would be natural to have an Error message. Please add two issues at Bitbucket for this :)

Comments

Thanks, and done (issues 347 & 348)