SIPG is misbehaving at the corners which separate a Dirichlet and a Neumann BC .

[Abbas]

Hello, 

Setting: 

1) Unit square domain, 
2) Simple scalar Laplacian with no source term.  
3) 0 and 1 Dirichlet on opposite ends of the square and zero Neumann on the other two edges. 

I am trying to solve the above problem using SIPG by adopting step 74. The expected solution should be perfectly linear with a constant gradient as the problem reduces to a pseudo 1-D problem. 

That is the case if I use CG, but with DG I do not get a perfectly linear solution and that is causing me problems.  

 
I have attached below a single image file that includes: 

1) The solution.

2) The solution gradient with weird behaviour at the corners 

3) Line plot of the solution across the lower wall. 

and the code I am using which is based on step 74. 

Things I tried/looked into: 

- I tried to adapt the code to CG and I get the expected results. I can't straight up impose BCs strongly with DG here. 

- I am sure this isn't a Paraview issue because I tried integrating across the boundary within deal and I am getting mesh dependent results. 
- The weak form is correct. The Neuman BC shows up as a boundary integral in the rhs which has no contribution. (Neumann zero at these walls).  
- I looked into the distributed LDG method within the dealii gallary and adopted it to these BCs and I didn't have this problem. 

The only thing I haven't tried yet is to do a simple problem by hand and compare the dealii and my outputs. Before doing that, any ideas on how I can debug this? Is this something associated with imposing BCs weakly with SIPG and will happen no matter what? 

best,

Abbas 

   

[Luca Heltai]

You have an issue in your boundary worker.


The boundary worker is called passing a cell, a face_no, scratch, and copy. Yet, in your code, you loop again on the faces of the cell:


for (const auto &face : cell->face_iterators())
{
if ( ((face->boundary_id() == 0) || (face->boundary_id() == 1)))


this should become

if((cell->face(face_no)->boundary_id() == 0) || (cell->face(face_no)->boundary_id() == 1) {

…

}

otherwise, you are going to assemble the face terms twice on the four corners.

Best,
Luca.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/df72814c-0cdc-4215-819c-46f5a8d05e84n%40googlegroups.com.
> <Screenshot from 2023-01-11 12-39-11.png><SIPG.cc>

[Abbas]

Well, this is embarrassing.
Thank you very much   

[Luca Heltai]

I didn’t mean to sound too harsh, sorry for that. I could spot the mistake right away because it is a mistake I do all the time.

No embarrassment is needed, and no question is ever a stupid question. Don’t worry! I’ve been in your shoes so many times I cannot count them anymore…

:)

Luca.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/35d379b6-4297-457e-90d3-14173b5c234dn%40googlegroups.com.

[Abbas]

I don't think you were harsh nor did take it that way. That was very considerate of you nonetheless. 

It was a light sleezy comment on my end. I should be more careful. 
Thanks again.