Problem with the use of FE_Nedelec with a mesh file generated from Gmsh

[Jianan Zhang]

Hi, all,

I am using dealii to solve the curl-curl equation:

curl(mu^(-1)curl(E)) + (-omega^2*epsilon+j*omega*sigma)*E=0,

with boundary conditions: n x E = n X F ,                          (Dirichlet)

                                    or  n X (curl(E)) = n X (curl(F))       (Neumann)

where F is a known exact solution.

Following Ross Kynch's code, I implemented my own exact solution and the assembly of system of equations my self. The exact solution (including real and imaginary parts) is as follows:

E_re = [0, sin(pi*x)/a*cos(k_z10_re*z)*exp(k_z10_im*z), 0],

E_im = [0, -sin(pi*x)/a*sin(k_z10_re*z)*exp(k_z10_im*z), 0]

The exact solution for Ross Kynch's case (take p=[0 1 0] and k=[0 0 0.1]) is given by:

E_re = [0, cos(0.1*z), 0],

E_im = [0, sin(0.1*z), 0] 

The domain I want to solve this problem is a 3D hyperrectangle  with width = 19.05e-3, height = 9.524e-3, and length = 15.24e-3. I've tried to solve the above equation with these two exact solutions. When I used the build-in mesh generator in dealii (i.e., GridGenerator::hyper_rectangle (triangulation, p1, p2);) to generate the mesh. The problem can be solved correctly. The results are given below:

my exact solution, built-in mesh generator:

Ross' exact solution, built-in mesh generator:

Cycle 0:Number of degrees of freedom: 3888

 Done

cycle cells dofs    L2 Error    H(curl) Error  

    0    64 3888 9.73783895e-12 1.03298196e-11 

When I use the mesh generated from Gmsh software (after conversion into hexes), if I give Ross' exact solution, the problem can still be solved correctly, as shown bellow:

Cycle 0:number of cells:244

Number of degrees of freedom: 13176

 Done

cycle cells dofs     L2 Error    H(curl) Error  

    0   244 13176 1.36020912e-07 8.58355818e-05

However, when I switch to my exact solution, the results are totally wrong. The electric field obtained has totally wrong direction (seems to be arbitrary but the correct one should be along y-axis). The results are as follows

   

I've spent like ten days finding the reason, and I am almost 100% sure that my exact solution is given correctly and the assembly is correct. The only thing I changed is the exact solution, and the difference between his and mine is that his is only z-component dependent while mine depend on both x and z. 

I post my codes and the vtk files for anyone who are interested. I am having  a hard time on debugging this and any suggestions and help is much appreciated. 

Thanks in advance.

Jianan Zhang

[Jianan Zhang]

Just forget to attach the mesh file generated using Gmsh.

Jianan

在 2017年12月20日星期三 UTC-5上午12:56:11，Jianan Zhang写道：

[Bruno Turcksin]

Jianan,

On Wednesday, December 20, 2017 at 12:56:11 AM UTC-5, Jianan Zhang wrote:

I've spent like ten days finding the reason, and I am almost 100% sure that my exact solution is given correctly and the assembly is correct. The only thing I changed is the exact solution, and the difference between his and mine is that his is only z-component dependent while mine depend on both x and z. 

I post my codes and the vtk files for anyone who are interested. I am having  a hard time on debugging this and any suggestions and help is much appreciated.

You cannot expect us to go through your code and debug it for you. We all have our own codes to debug... Now you say that the only difference between Ross' problem and yours is the dependence in both x and z. Now you have to determine if the problem is in the x component or in the x-z coupling. What happens if the solution depends only on x? only on y? What about y-z instead of x-z?

Best,

Bruno

[Jianan Zhang]

Hi, Bruno,

Thanks for the reply. Of course, I am not expecting someone to go through my code and debug it for me. I post the codes just in case that someone is interested or someone can take it as a reference when they try to use FE_Nedelec.

What I am really expecting is that someone who has some experience of using FE_Nedelec to solve EM problems can give me some general advice. For example, can we use FE_Nedelec  on a mesh imported from other mesh generators instead of the one generated from dealii itself? Or does dealii of the current version support this kind of mesh (obtained from Gmsh software and after transformed into hexahedrals) well? and so on and so forth.

If someone can give me this kind of hint it will help a lot. Because the actual problem may be due to the not good support of non-rectangular meshes when using FE_Nedelec, and in that case my debugging is like a wast of time.

Anyway, your suggestions make sense, I will try exact solutions of different forms and see when it works and when it does not work. Meanwhile, I will search previous works to see if anyone has met similar problems as me.

Best

Jianan

在 2017年12月20日星期三 UTC-5上午9:34:23，Bruno Turcksin写道：

[Bruno Turcksin]

Jianan,

2017-12-20 14:54 GMT-05:00 Jianan Zhang <adamzh...@gmail.com>:

What I am really expecting is that someone who has some experience of using FE_Nedelec to solve EM problems can give me some general advice. For example, can we use FE_Nedelec  on a mesh imported from other mesh generators instead of the one generated from dealii itself? Or does dealii of the current version support this kind of mesh (obtained from Gmsh software and after transformed into hexahedrals) well? and so on and so forth.

If someone can give me this kind of hint it will help a lot. Because the actual problem may be due to the not good support of non-rectangular meshes when using FE_Nedelec, and in that case my debugging is like a wast of time.

I don't follow closely what is happening with FE_Nedelec but it is probably the worst finite element to work with in deal. Over the years, we have found some subtle bugs in the implementation. My understanding is that the current implementation will be replaced by this one https://github.com/dealii/dealii/pull/2240 which should be a lot more robust. If you suspect that the problem is in deal, you could try this implementation.

Best,

Bruno

[Jianan Zhang]

Hi, Bruno,

Thanks for the prompt reply. I will try the implementation as you suggested, hope I can make it work with it. Otherwise, I will just go with brick elements, although the EM structure I can model may be limited.

Best

Jianan

在 2017年12月20日星期三 UTC-5下午3:26:21，Bruno Turcksin写道：

[Bruno Turcksin]

Jianan,

2017-12-20 17:21 GMT-05:00 Jianan Zhang <adamzh...@gmail.com>:

Thanks for the prompt reply. I will try the implementation as you suggested, hope I can make it work with it. Otherwise, I will just go with brick elements, although the EM structure I can model may be limited.

According to the comments in that PR it works on non-standard meshes. If your code doesn't work, I would contact the people working on the PR to make sure that your mesh works with the new implementation but don't exclude too quickly that the bug is your code ;-)

Best,

Bruno

[Jianan Zhang]

Bruno,

Okay, I will do more experiments and checks on my code and Ross' code to double confirm that the bug is not my code. Once I confirm this, it will be great if you can contact the people working on the PR to see whether my mesh can work with the new implementation. 

I'll contact you again in the near future!

Best

Jianan

在 2017年12月20日星期三 UTC-5下午5:35:11，Bruno Turcksin写道：

[Jianan Zhang]

Hi, Bruno, 

Hope you had a good holiday. I have done more checks on my code these days and found that the problem truly lies in the use of FE_Nedelec for non-rectangular meshes. 

To illustrate, I tried to solve the problem with Ross's exact solution (but with a larger value of omega) on the mesh generated by Gmsh, Then the results obtained are obviously wrong. But when I try the same exact solution on the mesh created by dealii, the results are correct. This demonstrates that dealii of the current version truly has not so good support for non-rectangular meshes.

So, can you help me contact the people working on the PR to try this out to see if my code works with the new implementation?

Many thanks

Jianan 

在 2017年12月20日星期三 UTC-5下午5:35:11，Bruno Turcksin写道：

Jianan,

[Bruno Turcksin]

Jianan

2017-12-28 10:40 GMT-05:00 Jianan Zhang <adamzh...@gmail.com>:

So, can you help me contact the people working on the PR to try this out to see if my code works with the new implementation?

If you just want to try the new implementation you can take Ross' branch and see if it works. What would be useful though is to create a test similar to this one but with a failing geometry. This way we can make sure that 1) it works without having to compile all of your code 2) it keeps working in the future by adding the test in our test suite.

Best,

Bruno

[Wolfgang Bangerth]

On 12/28/2017 08:40 AM, Jianan Zhang wrote:
>
> To illustrate, I tried to solve the problem with Ross's exact solution (but
> with a larger value of omega) on the mesh generated by Gmsh, Then the results
> obtained are obviously wrong. But when I try the same exact solution on the
> mesh created by dealii, the results are correct. This demonstrates that dealii
> of the current version truly has not so good support for non-rectangular meshes.

It's not deal.II itself, and it's not about non-rectangular meshes either --
it's a bug somewhere in the implementation of the Nedelec element. We've had
similar reports over the years before where something goes wrong with solving
on meshes that have non-standard face orientations. But we've never been able
to find exactly where in the ~4,000 lines of code that makes up the Nedelec
element the problem lies :-(

As Bruno noted, a small testcase is always helpful. I don't think we can
promise that anyone with the skills to find the bug has the time and
inclination to do so (it's a *volunteer* project, after all). But having a
small testcase is always a first step if anyone ever wanted to try to find the
problem.

Best
W.

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

[Jianan Zhang]

Thanks for your suggestions. I think I will first try the new implementation to see if it works. Later I may create a small testcase which can be added in the test suite.

Jianan 

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/ao-bi_Nh7gc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jianan Zhang]

Dear professor Bangerth,

Thanks for the reply. As I said, I will first try the new implementation myself to see if it works. Provided that it solves my problem correctly, I may create a small testcase that can be added into deal's test suite. 

Hope I can make some contribution in this regard.

Best

Jianan

[Jianan Zhang]

Hi, Bruno, 

Just want to inform you and the dealii users who may be interested that my problem is solved perfectly by using Ross' new implementation FE_NedelecSZ. The results given by brick elements and hexes are quite similar. Thanks a lot for the help.

Best

Jianan 

在 2017年12月28日星期四 UTC-5上午10:58:43，Bruno Turcksin写道：

[Bruno Turcksin]

Jianan,

2018-01-24 21:23 GMT-05:00 Jianan Zhang <adamzh...@gmail.com>:

Just want to inform you and the dealii users who may be interested that my problem is solved perfectly by using Ross' new implementation FE_NedelecSZ. The results given by brick elements and hexes are quite similar. Thanks a lot for the help.

Glad it worked and thank you for letting us know.

Best,

Bruno