Error using project_boundary_values_div_conforming with Raviart_Thomas FE

[t...@gmx.net]

Dear User Group,

I'm currently playing around with deal.ii to solve Maxwell's Equations and combining Tutorial pieces

I try to solve an Eigenvalue Problem (similar to the scalar-valued step-36) in Electrodynamics, namely find the Eigenmodes of a Spherical waveguide. The differential equation you try to solve there is the same as in the example in this document about Nedelec's Finite Elements (nedelec.pdf), namely curl(curl(E)) = omega^2*E. Now, for 2 spacial dimensions using Nedelec Elements with two in-plane Vector components you only calculate TE (Transversal Electric) modes, but this works fine. If I want to calculate the TM (Transversal Magnetic) modes, I THINK (please feel free to lecture me) I should use the differential equation for B, namely curl(curl(B)) = omega^2*B, but for B the normal components are required to be continuous, so I should use Raviart_Thomas Elements instead of Nedelec and use the Boundary condition of vanishing normal components, i.e. project_boundary_values_div_conforming(...). This is just changing 2 lines in the Code, but gives the error below at runtime (The Code is attached. As Error happens quite early in the make_grid_and_dofs function and the code is very much tutorial-like, I hope it is readable).
Additional question: is there a possibilty to use 3d-Nedelec Elements on a 2D geometry, or should I use an FE_System with a Nedelec Element and a FE_Q? In this case, It seems I have to reformulate the curl operator, which I'm not overly keen on...

Thanks in advance for your precious time!

Felix Schwarz

--------------------------------------------------------
An error occurred in line <77> of file </build/deal.ii-MRHzH9/deal.ii-8.4.2/source/fe/fe_poly_tensor.cc> in function
    void dealii::{anonymous}::get_face_sign_change_rt(const cell_iterator&, unsigned int, std::vector<double>&)
The violated condition was:
    f * dofs_per_face + j < face_sign.size()
The name and call sequence of the exception was:
    ExcInternalError()
Additional Information:
This exception -- which is used in many places in the library -- usually indicates that some condition which the author of the code thought must be satisfied at a certain point in an algorithm, is not fulfilled. An example would be that the first part of an algorithm sorts elements of an array in ascending order, and a second part of the algorithm later encounters an an element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an exception since it indicates an error in deal.II, not in your own program. Try to come up with the smallest possible program that still demonstrates the error and contact the deal.II mailing lists with it to obtain help.

Stacktrace:
-----------
#0  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2:
#1  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: dealii::FE_PolyTensor<dealii::PolynomialsRaviartThomas<2>, 2, 2>::fill_fe_values(dealii::TriaIterator<dealii::CellAccessor<2, 2> > const&, dealii::CellSimilarity::Similarity, dealii::Quadrature<2> const&, dealii::Mapping<2, 2> const&, dealii::Mapping<2, 2>::InternalDataBase const&, dealii::internal::FEValues::MappingRelatedData<2, 2> const&, dealii::FiniteElement<2, 2>::InternalDataBase const&, dealii::internal::FEValues::FiniteElementRelatedData<2, 2>&) const
#2  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: dealii::FEValues<2, 2>::do_reinit()
#3  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: void dealii::FEValues<2, 2>::reinit<dealii::DoFHandler, false>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > const&)
#4  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: void dealii::hp::FEValues<2, 2>::reinit<dealii::DoFHandler<2, 2>, false>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> >, unsigned int, unsigned int, unsigned int)
#5  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: void dealii::VectorTools::project_boundary_values_div_conforming<2>(dealii::DoFHandler<2, 2> const&, unsigned int, dealii::Function<2, double> const&, unsigned char, dealii::ConstraintMatrix&, dealii::Mapping<2, 2> const&)
#6  ./Microwave_Resonator: MR_fesc::EigenvalueProblem<2>::make_grid_and_dofs()
#7  ./Microwave_Resonator: MR_fesc::EigenvalueProblem<2>::run()
#8  ./Microwave_Resonator: main

[Wolfgang Bangerth]

Felix,

> I try to solve an Eigenvalue Problem (similar to the scalar-valued step-36) in
> Electrodynamics, namely find the Eigenmodes of a Spherical waveguide. The
> differential equation you try to solve there is the same as in the example in
> this document about Nedelec's Finite Elements (nedelec.pdf

> <http://www.math.chalmers.se/~stig/underv/doktorandkurs/FEM/200607/nedelec.pdf>),

> namely curl(curl(E)) = omega^2*E. Now, for 2 spacial dimensions using Nedelec
> Elements with two in-plane Vector components you only calculate TE
> (Transversal Electric) modes, but this works fine. If I want to calculate the
> TM (Transversal Magnetic) modes, I THINK (please feel free to lecture me) I
> should use the differential equation for B, namely curl(curl(B)) = omega^2*B,
> but for B the normal components are required to be continuous, so I should use
> Raviart_Thomas Elements instead of Nedelec and use the Boundary condition of
> vanishing normal components, i.e. project_boundary_values_div_conforming(...).

This can't be correct. The solutions of the curl-curl equation are in H_curl
and have continuous tangential component but possibly dicontinuous normal
components. You will not be able to impose normal boundary conditions.

> This is just changing 2 lines in the Code, but gives the error below at
> runtime (The Code is attached. As Error happens quite early in the
> make_grid_and_dofs function and the code is very much tutorial-like, I hope it
> is readable).

This is interesting, but independent of the issue above. Can you try to come
up with as small a program that shows this bug? If so, can you open a bug
report at
https://github.com/dealii/dealii/issues
for it?

> Additional question: is there a possibilty to use 3d-Nedelec Elements on a 2D
> geometry, or should I use an FE_System with a Nedelec Element and a FE_Q? In
> this case, It seems I have to reformulate the curl operator, which I'm not
> overly keen on...

I'm not sure I understand what you mean by "use 3d-Nedelec elements on a 2d
geometry". Can you elaborate?

Cheers
W.

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

[t...@gmx.net]

Wolfgang,

Thank you for your time.

>I THINK (please feel free to lecture me) I
> should use the differential equation for B, namely curl(curl(B)) = omega^2*B,
> but for B the normal components are required to be continuous, so I should use
> Raviart_Thomas Elements instead of Nedelec and use the Boundary condition of
> vanishing normal components, i.e. project_boundary_values_div_conforming(...).

This can't be correct. The solutions of the curl-curl equation are in H_curl
and have continuous tangential component but possibly dicontinuous normal
components. You will not be able to impose normal boundary conditions.

 

Ok, thank you for that comment, I'd better do some more thinking on that then.

> This is just changing 2 lines in the Code, but gives the error below at
> runtime (The Code is attached. As Error happens quite early in the
> make_grid_and_dofs function and the code is very much tutorial-like, I hope it
> is readable).

This is interesting, but independent of the issue above. Can you try to come
up with as small a program that shows this bug? If so, can you open a bug
report at
   https://github.com/dealii/dealii/issues
for it?

 

Done, this is issue 4092 now.

> Additional question: is there a possibilty to use 3d-Nedelec Elements on a 2D
> geometry, or should I use an FE_System with a Nedelec Element and a FE_Q? In
> this case, It seems I have to reformulate the curl operator, which I'm not
> overly keen on...

I'm not sure I understand what you mean by "use 3d-Nedelec elements on a 2d
geometry". Can you elaborate?

Can I calculate a 3d-Vector Field that only depends on two spacial Koordinates, in this case (E_x(x,y),E_y(x,y),E_z(x,y))? That, I think, means E is a function E: R^2->R^3. Even more general: in the Documentation about the VectorTools::project_boundary_values_curl_conforming_l2 function, there is something written about a six-component Vektor-Field, combining E and B, in 3d space. I did not find any hints on how to actually write that other than with FE_System.

Regards,

Felix

[Wolfgang Bangerth]

> This is interesting, but independent of the issue above. Can you try to come
> up with as small a program that shows this bug? If so, can you open a bug
> report at
> https://github.com/dealii/dealii/issues
> <https://github.com/dealii/dealii/issues>
> for it?
>
>
> Done, this is issue 4092 now.

Great, thanks. I can reproduce this, but I don't know when I (or anyone else)
may be able to get around to fixing this :-(

Best