using VectorTools::project_boundary_values_curl_conforming_l2 with FEsystem fe(FE_Nedelec<3>(0), 2)

[Abbas Ballout]

Some context: 
I am trying to model a waveguide which would require me to solve the double curl equation curlxcurlxE-cE=0 where c is some constant(real or imaginary) subject to nxE=nxf on different parts of the boundary where f can be real or imaginary also. 

The problem requires the use of Nedelec elements but you already know that. Since the problem is complex I am breaking it down into two with FEsystem. I am trying to approach this by solving this simple system before complicating things:
curlxcurlxE_r-cE_r=0
nxE_r=1   on an "inlet"
nxE_r=0   everywhere else 
and 
curlxcurlxE_i-cE_i=0
nxE_i=1   on an "inlet"
nxE_i=0   everywhere else 

These two systems are uncoupled and exactly the same. The issue has been that my BCs aren't being applied for the second part of the equation. I am using VectorTools::project_boundary_values_curl_conforming_l2(). The first part of the equation gives a solution (inaccurate but that's not my issue here) with the BCs applied. 

I print out the sparsity pattern and there seems to be a difference between applying BCs and not applying them which makes me conclude that project_boundary_values_curl_conforming_l2 is effecting the matrix. 
I printed out the global right hand side for a smaller problem and the second half of it, which is the half that is related to the second problem, is just zeros so my guess is that constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs) isn't applying the BC in the cell_rhs when FENedelec systems are involved. 

I have attached my code here.
Thanks  

[Wolfgang Bangerth]

On 5/2/21 10:41 AM, Abbas Ballout wrote:
>
> I print out the sparsity pattern and there seems to be a difference between
> applying BCs and not applying them which makes me conclude that
> project_boundary_values_curl_conforming_l2 is effecting the matrix.
> I printed out the global right hand side for a smaller problem and the second
> half of it, which is the half that is related to the second problem, is just
> zeros so my guess is that constraints.distribute_local_to_global(cell_matrix,
> cell_rhs, local_dof_indices, system_matrix, system_rhs) isn't applying the BC
> in the cell_rhs when FENedelec systems are involved.

These functions should not affect the sparsity pattern, but only the matrix.
Indeed, they don't consider the matrix themselves, but only build the constraints.

It is entirely possible that the function you are calling has a bug, but
before we conclude this, I would suggest that you output the constraints each
of the four calls you have in your code generates. After all, this is what
this function produces, and you should verify that these are correct -- maybe
using a mesh that only consists of one cell for which you can verify correctness.

Best
W.

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

[Abbas Ballout]

What I meant by a sparsity pattern change was that the system becomes less dense when the BCs are applied Sorry about that. 
I will elaborate on my previous mail here because I don't believe that I actually did a good job at it. 
So I did 4 different simulations for the same problem as above ( FEsystem with 2 Nedelecs ) for a 4 cell cubic 3D geometry. 

1. The first simulation has no BCs applied
2. The second simulation has BCs applied only for the first part of the system (only the first equation)
3. The third simulation has BCs applied only for the second part of the system (only the third equation)
4. The third simulation has BCs applied to both equations.

and I print out the sparsity pattern and the system right hand side (I attached both down below) and here's what the problem seems to be:

The right hand side of the second equation isn't being updated when I apply BCs even though the system matrix is. 

1) Hints on what I might be missing in my code? 
2) Any suggested workarounds in case that there is an actual bug?

[Daniel Arndt]

Abbas,

There was a seemingly related issue fixed in deal.II 9.2, see https://github.com/dealii/dealii/issues/8974 and https://github.com/dealii/dealii/pull/9000.

Are you using the latest release and testing in Debug mode?

Do you have a minimal example reproducing the issue to share?

Best,

Daniel

--
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/a599c2d5-8bf2-4054-9d10-5555590829f1n%40googlegroups.com.

[Wolfgang Bangerth]

On 5/3/21 7:43 AM, Abbas Ballout wrote:
> *The right hand side of the second equation isn't being updated when I apply

> BCs even though the system matrix is.

> *

> 1) Hints on what I might be missing in my code?
> 2) Any suggested workarounds in case that there is an actual bug?

That's too far downstream. You have a suspicion that a function isn't working
as expected (which is totally plausible), so let's check that right after that
function is called. Dealing with the right hand side happens much later in
your program, and many other things could have gone wrong in the meantime. I
would just print the constraints right after they were computed and check that
first.

Best
W>

[Abbas Ballout]

I am running my code in Debug using the latest 9.3 release. 

I have attached a minimal code and what it does it that it prints the constraints.

• If BCs are applied for both parts of the system ,my output for the constraint is:

1.     0 = 0
2.     1 = 0
3.     2 = 0
4.     3 = 0
5.     4 = 1.00000
6.     5 = 0
7.     6 = 1.00000
8.     7 = 0

• If BCs are applied for the first part of the system, my output for the constraint is:

1.     0 = 0
2.     2 = 0
3.     4 = 1.00000
4.     6 = 1.00000

• If BCs are applied for he first partof the system my output for the constraint is:

1.     1 = 0
2.     3 = 0
3.     5 = 0
4.     7 = 0

[Wolfgang Bangerth]

On 5/4/21 6:15 PM, Abbas Ballout wrote:
> I am running my code in Debug using the latest 9.3 release.
> I have attached a minimal code and what it does it that it prints the constraints.
>

> * If BCs are applied for both parts of the system ,my output for the
> constraint is:
>
> 1.     0 = 0
> 2.     1 = 0
> 3.     2 = 0
> 4.     3 = 0
> 5.     4 = 1.00000
> 6.     5 = 0
> 7.     6 = 1.00000
> 8.     7 = 0
>
> * If BCs are applied for the first part of the system, my output for the
> constraint is:
>
> 1.     0 = 0
> 2.     2 = 0
> 3.     4 = 1.00000
> 4.     6 = 1.00000
>
> * If BCs are applied for he first partof the system my output for the
> constraint is:
>
> 1.     1 = 0
> 2.     3 = 0
> 3.     5 = 0
> 4.     7 = 0
Abbas,
excellent little test case -- thanks for shrinking things to a size that makes
it easy for me to understand what is happening! I'm going to add a variation
of your testcase to the test suite!

The problem is that you have 6 finite element components, but your boundary
value function only provides 3 components. The function
VectorTools::project_boundary_values_curl_conforming_l2 was missing a check
for this to tell you about this requirement.

In any case, you can make things work by using this instead:
```
class BC_class : public Function<3>
{
public:
BC_class()
: Function<3>(6)
{}

void
vector_value_list(const std::vector<Point<3>> &points,
std::vector<Vector<double>> &values) const override;
};

void
BC_class::vector_value_list(const std::vector<Point<3>> &points,
std::vector<Vector<double>> &values) const
{
for (unsigned int i = 0; i < points.size(); ++i)
{
values[i][0] = 1.0;
values[i][1] = 0.0;
values[i][2] = 0.0;
values[i][3] = 1.0;
values[i][4] = 0.0;
values[i][5] = 0.0;
}
}
```

Best
W.