Applying PointSource fails with 1 element in 1D

[Jack Hale]

A bit of an edge case, but applying a PointLoad to the end node of the 1D mesh fails. I am running Dolfin 1.0 so I guess this may have been fixed since.

from dolfin import *

mesh = UnitInterval(1)

V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

end_point = Point(1.0)
f = PointSource(V, end_point)

a = inner(u,v)*dx
L = Constant(0.0)*v*dx

A_matrix = assemble(a)
b_vector = assemble(L)

f.apply(b_vector)

with the following error:

Traceback (most recent call last):
  File "bug_pointsource.py", line 18, in <module>
    f.apply(b_vector)
StandardError: CGAL ERROR: assertion violation!
Expr: m_primitives.size() > 1
File: /usr/include/CGAL/AABB_tree.h
Line: 302

[Jack Hale (jack-hale)]

I actually just noticed this bug is very similar (possibly the same) as the one in Prof. Arnold's post:

https://bugs.launchpad.net/dolfin/+bug/1040312

I will compile the latest version from bzr and see if it is fixed.

[Jack Hale (jack-hale)]

OK, I compiled the latest version using dorsal with STABLE_BUILD=false and I am getting the same issue.

The fix commited by Johan Hake:

http://bazaar.launchpad.net/~dolfin-core/dolfin/trunk/revision/6859

for Prof. Arnold's bug:

https://bugs.launchpad.net/dolfin/+bug/1040312

doesn't seem to be helping here.

[Johan Hake (johan-hake)]

The bug _is_ related to:

  https://bugs.launchpad.net/dolfin/+bug/1040312

The fix applied to the former bug was to specific to help you. One can argue that the former bug fix was also implemented at the wrong place. It should have been implemented either in the Mesh class or in the IntersectionOperator, where it would have helped you...

[Jack Hale (jack-hale)]

OK, so I can see currently that the PointSource::apply() method calls the Mesh::intersected_cell() method which in turn will call the appropriate method in the appropriate IntersectionOperator class which is dependent on the primitive type of the mesh.

So your suggestion is to move the logic you submitted in the patch to the Mesh::intersected_cell(const Point& point) method so that it will be used consistently across DOLFIN.

[Johan Hake (johan-hake)]

Yes. Should be easy enough. But the problem is that there are many
Mesh::intersected_foo, (and some Mesh::closest_foo) which need to be
covered by the logic...

Any patches or merge requests with added unit tests are most welcomed ;)!

Johan

On 08/30/2012 07:02 PM, Jack Hale wrote:
> OK, so I can see currently that the PointSource::apply() method calls
> the Mesh::intersected_cell() method which in turn will call the
> appropriate method in the appropriate IntersectionOperator class which
> is dependent on the primitive type of the mesh.
>
> So your suggestion is to move the logic you submitted in the patch to
> the Mesh::intersected_cell(const Point& point) method so that it will be
> used consistently across DOLFIN.
>

[Jack Hale (jack-hale)]

Yep I can give this a go... I haven't signed the copyright release form yet though.

[Jack Hale (jack-hale)]

OK, I've put the new fix in my work branch on launchpad. I have run it against both my example and Prof. Arnold's example and it fixes both CGAL errors.

I have put the copyright release form in the post today, I consent to the terms.

[Johan Hake (johan-hake)]

Thanks!

I made some small changes and added similar fixes for some more methods.

Johan

On 09/03/2012 01:22 PM, Jack Hale wrote:
> OK, I've put the new fix in my work branch on launchpad. I have run it
> against both my example and Prof. Arnold's example and it fixes both
> CGAL errors.
>
> I have put the copyright release form in the post today, I consent to
> the terms.
>

[Johannes Ring (johannr)]

This was fixed in 1.1.0.