Neumann boundary conditions and vanishing surface integrals, 3 questions

variational formulation
1

Hello FEniCS community!
I have several questions regarding the implementation of the Neumann boundary conditions in a very specific case.

I am dealing with a mixed function space that models, say temperature and voltage (2 scalar fields). The domain is a 3D box, where I set a Dirichlet boundary condition on a subsurface of a surface of that box. I set a Neumann boundary conditions on a part of another surface of that box. This is enough to ensure that the solution exists and is unique through that box. The full weak form is extremely long and complicated, and there are several surface terms that, according to the FEniCS tutorial, vanish at the surfaces where Dirichlet boundary conditions are applied.

My 1st question is whether it hurts if I do consider those surface integrals arising when integrating by parts second order spatial derivatives, knowing that they should vanish on the (sub)surfaces where Dirichlet boundary conditions are set. Or should I really avoid them (in an ad-hoc way as the tutorial does) as if they were equal 0. My experience is that it hurts a lot if I do not throw out those surface integrals, but that seems so strange to me that this leads me to ask this question here. Because after all, if they really vanish, FEniCS/dolfin should evaluate them as close to being 0, it shouldn’t hurt that much. But that’s definitely not what I get.

My second question is whether any kind of surface integral should vanish on those surfaces? I guess not. For example if I set the temperature on such a subsurface to be fixed as a Dirichlet boundary condition, then only surface integrals dealing with the function space of temperature should vanish there, but not the one dealing with the voltage function space. Is that correct?
What I mean by that, is that if

boundary_parts = MeshFunction("size_t", mesh, some_file)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
u, v = TestFunctions(ME)           
TV = Function(ME)
Temperature, Voltage = split(TV)
temperature_value = 100.0
bc = [DirichletBC(ME.sub(0), temperature_value, boundary_parts, fixed_temperature_surface)]

then surface integrals having the term “u” should vanish on the “fixed_temperature_surface” surface, but not the surface integrals having the term “v”. Is that correct? In my experience, this is correct, but I just want to make sure.

My third question is, what happens with the surface integrals on the rest of the box? I.e. where I haven’t specified any Neumann nor Dirichlet boundary conditions? I’ve read that if nothing is specified, Neumann boundary conditions are set by default. This would mean that I should consider almost all the surface of the box with the surface integrals. In fact the whole surface of the box for the surface integrals that deal with “v”. That would make sense if I were to follow the tutorial to the letter. However from my experience, this leads to non nonsensical solutions. What leads to sensical solutions is to consider only the surface(s) where I explicitly set the Neumann boundary conditions. Not the implicit ones. I would really, really appreciate a confirmation regarding this remark.

Thank you!

2

FEniCS does not know that you have Dirichlet BCs on a surface before you apply the BCs, so what happens is:

  1. Assemble all rows of the system matrix, also those related to DOFs on the Dirichlet surface
  2. Apply Dirichlet BCs by “removing” rows in the matrix whos DOFs are on the Dirichlet boundary and replace them by rows with “1” on the diagonal, zero everywhere else and the Dirichlet value in the same location on the right hand side. Dolfin may optionally “remove” also the column containing these DOFs (maintaining symmetry)

So, the integrals do not “vanish”, they may in fact give very large contributions, it is just that the row in the matrix is repurposed to set strong BCs afterwards and then the contribution is forced to zero.

PS: you can apply Dirichlet BCs weakly (Nitsche’s method), and you may get ever so slightly better results since you make the linear solver try to fulfill both the PDE and the BCs at the boundary. Strong BCs fulfill only the BCs at the boundary.

PPS: Neumann BCs “by default” is a common misconception. Often Neumann BCs with derivative equal to zero result in surface integrals that are identically zero, so not including them or including them are the same. But, if the Neumann surface integrals do not evaluate to zero then it is obviously wrong to drop them.

2 Likes