Dirichlet boundary conditions implementation bug

[Vojtech Kulvait]

I noticed very strange things regarding Dirichlet boundary conditions when debugging my code. I defined boundary conditions but it seems that it wont apply them on two particular points which is a problem also when obtaining solution.

I have a bit special domain and I am enclosing minimal python source code to reproduce the bug.

After running reproduceBug("xmlfiledir") i suppose to see u2.vector()==u1.vector() but I do see u1.vector()[916]==1 u2.vector()[916]==0

[Jan Blechta]

This report is ill-posed. Nobody can follow your 164 lines of confusing code to infer that u2.vector()==u1.vector() should hold. Moreover, there is no guarantee of some DOF ordering generally.

[Vojtech Kulvait]

-

[Vojtech Kulvait]

Finally it is a bug! :) The problem is that DirichletBC class computes boundary conditions on facet level not DOF level. So if there are some nodal values with one boundary condition and other nodal values with different boundary condition sharing same facet the behavior is following: precedence takes the boundary condition, where the facet lies fully in the BC specification and is applied last to the whole facet. This may mean there are some nodal points without BC being even applied. I do think boundary conditions should be applied to all dofs nodal values individually. Function SubDomain::apply_markers should not contain all_points_inside and probably the whole DirichletBC::apply should only tabulate all nodal values and then test each of them against SubDomain::inside. I know there can be some speed cost but at least it will be correct.

[Vojtech Kulvait]

BUG HERE

[Martin Sandve Alnæs]

Good that you understand what happens, but this is expected behaviour and not a bug. Read the docstring of DirichletBC, in particular about the 'method' argument.

[Vojtech Kulvait]

So I am sorry. You are right:) So finally was the problem on my side. Sometimes I have hard time to go through. Thank you.

[Vojtech Kulvait]

You can reset to invalid I could not do that myself and did not catch your action.

[Vojtech Kulvait]

Well the problem is that even if I change the method "pointwise". The statement // Check if the coordinates are part of the sub domain (calls // user-defined 'inside' function) Array<double> x(gdim, &data.coordinates[i][0]); if (!_user_sub_domain->inside(x, false)) continue; It is on 1050 row in https://bitbucket.org/fenics-project/dolfin/src/ed96da8013c30b265fc81a785243e4c5c7e05c1b/dolfin/fem/DirichletBC.cpp?at=master

ends the loop because my bc returns (on_boundary && dolfin::near(x[0], 1.0)) or something similar. So the result is that the boundary condition is not applied anyway.

The only way is to check in subdomain implementation wether it lies on boundary and do not rely on on_boundary value. If this is a feature I am leavin this thread forever alone...

[Jan Blechta]

Yes, on_boundary does not work for method='pointwise'. The documentation does not reflect it.

[Prof Garth Wells]

Fix Issue #309.

→ <<cset 37e1a0bfa54b>>

[Shankha Nag]

Hi,

I have the similar problem. I am using method 'pointwise' and not using on_boundary. But still some of the nodes are not set BC.

I have to set different BC to 72 nodes, out of which 69 are set properly and 3 are not set.

Any expect advice is highly welcome.

Sincerely, Shankha