MappingQ, and minimal test for reproducing Benchmark

[Juan Carlos Araujo Cabarcas]

Dear all,

I'm am trying to reproduce with my implementation, the results in the Photonic Crystal computations performed in [1]. Here, the author uses a grid with an inner disk with radius R=0.475, and for FEM it is used the software Concepts [2] that implements curvilinear elements denoted Blending technique, or transfinite interpolation in quadrilaterals [3], [4,p. 144]. Exponential convergence was obtained for the p-FEM version.

By using the deal.II SphericalManifold, I can only get with the attached code up to R<0.46, before I get the exception:

"The image of the mapping applied to cell with center [0.446317 -0.206317] is distorted. The cell geometry or the mapping are invalid, giving a non-positive volume fraction of -0.000129188 in quadrature point 54."

I am attaching a minimal test case (modified from step-10) that reproduces the errors:
    mesh_test.cc
and the code for the same mesh used in [1], adapted to deal.II standards in:
    unit_cell.h

The question is then if it is possible to get this example running for R=0.475, in order to reproduce the results in [1]. Is there something I am doing wrong that can be improved? Can this behaviour be explained by round off's?

Furthermore, is anyone implementing transfinite interpolation [3,4] in deal.II?

Thanks in advance, any help is greatly appreciated!

References:
[1] Engström, C.,Spectral approximation of quadratic operator polynomials arising in photonic band structure calculations, Numer. Math. (2014) 126: 413.
[2] Concepts web page: http://wiki.math.ethz.ch/Concepts/WebRoot?redirectedfrom=Concepts.WebHome
[3] W. J. Gordon, C. A. Hall, Construction of curvilinear coordinate systems and application to mesh generation Int. J. Num. Meth. Eng., Vol. 7 (1973), pp. 461-477
[4] Pavel Solin, Karel Segeth, Ivo Dolezel, Higher-Order Finite Element Methods, Studies in Advanced Mathematics, CRC Press, 2003.

[Martin Kronbichler]

Dear Juan Carlos,

We have now finalized the TransfiniteInterpolationManifold implementation in deal.II, 

http://dealii.org/developer/doxygen/deal.II/classTransfiniteInterpolationManifold.html

If you are interested, please try it out and let us know. As its implementation is completely new and I have only applied it to two applications in my group, there might still be issues, so let us know.

By the way, if you are working on a 2D disk I recommend using PolarManifold rather than SphericalManifold. I have not yet looked in your code, but I recommend to try the transfinite interpolation manifold first.

Best,
Martin

[Juan Carlos Araujo Cabarcas]

Thanks a lot Martin,

I had been waiting for this type of solution in dealii for a long time!
I also read in another thread about using MappingQGeneric instead of MappingQ, so I updated my test code as below, and everything worked like a charm!

                    PolarManifold<dim> boundary;
                    TransfiniteInterpolationManifold<dim> inner_manifold;
                   
              Triangulation<dim> triangulation;
             
                    Geom_parameters gp;
                    unsigned int curved_faces_label=10, not_curved = 1;
                   
                    cell_rod ( triangulation, 1.0, R, gp, false );
                   
                    triangulation.set_manifold ( curved_faces_label, boundary ); // labels defined with the grid
             
              inner_manifold.initialize(triangulation);
              triangulation.set_manifold (not_curved, inner_manifold);
             
              const MappingQGeneric<dim> mapping (degree);
                    const QGauss<dim> quadrature(degree);
             
              // const FE_Q<dim>     dummy_fe (1);
                    const FE_Q<dim> dummy_fe (QGaussLobatto<1>(degree + 1));
              DoFHandler<dim> dof_handler (triangulation);

              FEValues<dim> fe_values (mapping, dummy_fe, quadrature,
                                       update_JxW_values);

Juan Carlos Araújo,
Umeå Universitet