- changed status to invalid
Dirichlet boundary conditions implementation bug
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
Comments (12)
-
-
reporter -
-
reporter 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.
-
reporter - changed status to open
BUG HERE
-
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.
-
- changed status to invalid
-
reporter - changed status to resolved
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.
-
reporter You can reset to invalid I could not do that myself and did not catch your action.
-
reporter - changed status to open
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...
-
Yes,
on_boundary
does not work formethod='pointwise'
. The documentation does not reflect it. -
- changed status to resolved
Fix Issue
#309.→ <<cset 37e1a0bfa54b>>
-
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
- Log in to comment
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.