Point integrals

I'd like to be able to use 'point' measures the same way as I can use 'boundaries' (ie: not have them limited to PointIntegralSolver).

Issue from this question: http://fenicsproject.org/qa/4097/point-integrals?show=4104#a4104

eg:

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)

Comments (5)