A problem about cylindrical coordinates.

[RUIJIAN LU]

Hi everyone, 
I am a beginner to dealii, now I want to solve an equation on a cylindrical coordinate system. I create a triangulation by cylinder() and use CylindricalManifold for it, but it has some error. How can I do with it?

error:

Number dealii::Physics::VectorRelations::signed_angle(const dealii::Tensor<1, dim, Number>&, const dealii::Tensor<1, dim, Number>&, const dealii::Tensor<1, dim, Number>&) [with int spacedim = 3; Number = double]
The violated condition was:
    std::abs(axis * a) < 1.e-12 * a.norm() && std::abs(axis * b) < 1.e-12 * b.norm()
Additional information:
    The vectors are not perpendicular to the axial vector.

code:

    const Tensor<1, 3> axis({1.0, 0.0, 0.0});

    const Point<3> axial_point(0.0, 0.0, 0.0);

    const CylindricalManifold<3> cylinder(axis, axial_point);

   

    GridGenerator::cylinder(triangulation, 2.5, 2.5);

    for (auto cell : triangulation.active_cell_iterators())

        cell->set_all_manifold_ids(2);

    triangulation.set_manifold(2, cylinder);

[RUIJIAN LU]

Oh, the error occurs after calling triangulation.refine_global(2);
I

[Wolfgang Bangerth]

On 3/22/23 01:47, RUIJIAN LU wrote:
>
> I am a beginner to dealii, now I want to solve an equation on a cylindrical
> coordinate system. I create a triangulation by cylinder() and use
> CylindricalManifold for it, but it has some error. How can I do with it?
>
> error:
> Number dealii::Physics::VectorRelations::signed_angle(const dealii::Tensor<1,
> dim, Number>&, const dealii::Tensor<1, dim, Number>&, const dealii::Tensor<1,
> dim, Number>&) [with int spacedim = 3; Number = double]
> The violated condition was:
>     std::abs(axis * a) < 1.e-12 * a.norm() && std::abs(axis * b) < 1.e-12 *
> b.norm()
> Additional information:
>     The vectors are not perpendicular to the axial vector.

Can you show what the backtrace is that gets you to this place? It would be
important to understand what these vectors a, b, and axis are.

Best
W.

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

[Marc Fehling]

Hi!

The manifold should be attached only to the hull of the cylinder, but in your code snippet you attach it to all cells. This leads to problems when new vertices need to be created on the symmetry axis.

After a look in the GridGenerator::cylinder() documentation, you will notice that a CylindricalManifold object is already attached to those faces with a manifold ID 0. So you don't have to attach an additional one yourself!

Marc

[RUIJIAN LU]

Hi Marc!
Thanks, you are right. Therefore, when I just create a triangulation by cylinder(), have the coordinates of its vertices already been in cylindrical coordinate（r, phi, z） system？(Of course, the 'z' should be x.) If it is right, is fe_values.quadrature_point(q_index)(0) equal to ‘r’, fe_values.quadrature_point(q_index)(1)equal to 'phi' and fe_values.quadrature_point(q_index)(0) equal to 'x'? 
Best,

LU

[RUIJIAN LU]

Hi Dr. Bangerth,

Thanks for your message. 

I think the error is like Marc said. And if I call triangulation.refine_global() before set_manifold(), the error will go away. 

Best,
LU

[Wolfgang Bangerth]

On 3/23/23 10:32, RUIJIAN LU wrote:
>
> Thanks, you are right. Therefore, when I just create a triangulation by
> cylinder(), have the coordinates of its vertices already been in
> cylindrical coordinate（r, phi, z） system？(Of course, the 'z' should
> be x.) If it is right, is fe_values.quadrature_point(q_index)(0) equal
> to ‘r’, fe_values.quadrature_point(q_index)(1)equal to 'phi' and
> fe_values.quadrature_point(q_index)(0) equal to 'x'?

deal.II only ever computes in a Cartesian coordinate system. Let's call
it x1,x2,x3. If you would like to *interpret* x1 as the radial
direction, and x2 as the angle phi, and x3 as the vertical distance z,
then that is your *choice* when you write down the equation. But the
coordinates of the cylinder's vertices are stored in x1,x2,x3.

[RUIJIAN LU]

Hi Dr. Bangerth,

Thanks for your message. 

I understand it. So when computing the equation and I need it times "2*pi*r", I should call pull_back() to transform the coordinates of the cylinder's vertices to "(r, phi, z)" to get "r". Is it right?

Best,

LU

[Wolfgang Bangerth]

On 3/23/23 22:36, RUIJIAN LU wrote:
>
> I understand it. So when computing the equation and I need it times
> "2*pi*r", I should call pull_back() to transform the coordinates of the
> cylinder's vertices to "(r, phi, z)" to get "r". Is it right?

First, if you work in cylindrical coordinates, you probably have a 2d
domain in (r,z) space. It is probably rectangular.

deal.II only knows that you have a domain in x1,x2 space. If you ask a
cell for a vertex, for example, it will give you the coordinates in
x1,x2 space. For yourself, you will then say "the x1 coordinate
corresponds to my r coordinate" and "the x2 coordinate corresponds to
the z coordinate".

So where you need 2*pi*r in an integration formula, you will do
something like
2*pi*fe_values.quadrature_point(q)[0] // [0] == x1 coordinate == r