interpolate_boundary_values not adding to ConstraintMatrix

[mrjonm...@gmail.com]

Hi,

I'm working on a flux-pressure model using a Raviart-Thomas/DGQ FESystem. I'm using Dirichlet conditions on the RT (flux) elements on part of the domain and on the DGQ (pressure) elements on another part of the domain using the boundary_id's. I'm using VectorTools::interpolate_boundary_values to add the pressure BCs, but it doesn't appear to be adding rows to the ConstraintMatrix. 

The enclosed minimal code shows that interpolate_boundary_values is computing the BC value function, but it also shows that there is no change in the ConstraintMatrix after it is called. I can't figure out what logic in do_interpolate_boundary_values is keeping it from adding the appropriate entries to the ConstraintMatrix.

I would appreciate if anyone can find my error here or explain why I'm not getting the constraints I'm expecting. 

Thanks,

Jonathan

[Wolfgang Bangerth]

On 11/27/18 9:02 PM, mrjonm...@gmail.com wrote:
> I'm working on a flux-pressure model using a Raviart-Thomas/DGQ FESystem. I'm
> using Dirichlet conditions on the RT (flux) elements on part of the domain and
> on the DGQ (pressure) elements on another part of the domain using the
> boundary_id's. I'm using VectorTools::interpolate_boundary_values to add the
> pressure BCs, but it doesn't appear to be adding rows to the ConstraintMatrix.

That's correct -- and it is because the DG elements have no degrees of freedom
on faces, they are all located in the interior of the cell.

In other words, you can't impose boundary conditions strongly on DG elements.
You need to impose them in your variational formulation, for example using a
penalty term in the same way as you are presumably adding penalty terms in the
interior of the domain. (Unless, of course, you use a formulation that has no
derivatives on the pressure variable at all -- in which case the Dirichlet
values on the pressure should show up in the bilinear form to begin with;
isn't this how step-20 does it?)

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[mrjonm...@gmail.com]

Thanks for that explanation. I'll use the boundary integral like in step-20. 

When you say "the DG elements have no degrees of freedom on faces, they are all located in the interior of the cell", I understand how that's true for degree 0, but I'm guessing it's only in the "logical" sense for higher degree elements? I was just surprised that I still didn't get any constraints out of it for order 1 and 2 elements (even though I now know that's not the right thing to do.)

[Wolfgang Bangerth]

On 11/28/2018 10:48 AM, mrjonm...@gmail.com wrote:
>
> When you say "the DG elements have no degrees of freedom on faces, they
> are all located in the interior of the cell", I understand how that's
> true for degree 0, but I'm guessing it's only in the "logical" sense for
> higher degree elements?

Correct. We say that a degree of freedom is "logically" located on a
face if the shape function is continuous across that face in some sense.
That's because the face is shared between two cells, and so the DoF is
as well. But for DG elements, each DoF is uniquely associated with just
one cell, so it is not logically located on a face (or vertex, or edge).

It is true that DG elements are implemented as Lagrange interpolation
elements with interpolation points that are physically located on the
face. But logically they are not. (And we could have implemented these
elements in many other ways as well -- for example by having the
interpolation points located at Gauss points that really truly are in
the interior of the cell.)

> I was just surprised that I still didn't get any
> constraints out of it for order 1 and 2 elements (even though I now know
> that's not the right thing to do.)

I hope the explanation above makes clear why -- there are no degrees of
freedom that are logically located on the boundary; they are logically
all located in the interior, and consequently trying to compute boundary
values is not possible.