system_to_component_index() usage in face_worker for DG method

[Sylvain Mathonnière]

Hello everyone !

I have (another) a problem in my code. I am trying to mix step-12 and step-31 for my project.
I want to solve two coupled equations (named eq1 and eq2) that I am solving in parallel just like in step-31 but for one of them I am using the Discontinuous Galerkin method like in step-12 (eq1) and for the other one I use standard Galerkin (eq2). I also have two dof_handler like in step-31 one for each equation. It looks like this :

         FESystem<dim>        fe_eq1;         //FE element for equation 1
         FESystem<dim>        fe_eq2;         //FE element for equation 2
         DoFHandler<dim>      dof_handler_eq1;
         DoFHandler<dim>      dof_handler_eq2;

And the constructor of the class :

      template <int dim> DG_FEM<dim>::DG_FEM ()
        : fe_eq1 (FE_DGQ<dim>(1),40), fe_eq2 (FE_Q<dim>(1), 1), dof_handler_eq1             (triangulation), dof_handler_eq2 (triangulation) {}

For the eq1, I am using like in step-12 the MeshWorker::mesh_loop(...) to fill my matrix :


In the cell_worker function, in order to access the ith component of my shape function, I use the same technique as in step-8:

        ...
        
        const unsigned int n_dofs_per_cell = scratch_data.fe_values.get_fe().dofs_per_cell;

        for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
        {
            unsigned int component_i = fe_eq1.system_to_component_index(i).first;

            // and then I access the "component_i" of my shape function "i" at gauss point "point" using

            copy_data.cell_rhs(i) += .... fe_v.shape_value_component(i, point, component_i) ...
        }
        
        ...

This seems to work fine. Since I have 40 components and 4 nodes, I have n_dofs_per_cell = 160. That is confirmed by the debuggeur.

My problem arises when I am using the same technique in the face_worker.
        
        ...

        const unsigned int n_dofs        = fe_iv.n_current_interface_dofs();

        for (unsigned int i = 0; i < n_dofs; ++i)
        {
          unsigned int component_i = fe_eq1.system_to_component_index(i).first;
   
          copy_data_face.cell_matrix(i, j) += * fe_iv.shape_value((beta_dot_n > 0), j, qpoint, component_i)
        
        }
        
        ...

Here n_dofs = 320 and so system_to_component_index() does not like to have an argument i>=160 and throws an error.

How come n_dofs=320 in this case ? I was expecting 2 interfaces * 2 nodes * 40 components = 160 dofs, unless I am considering all the interfaces of the cell at onces in which case I have 2 interfaces * 4 nodes * 40 components = 320.

My main question is therefore:

How can I access the ith component of my shape function in face_worker like in the cell_worker if I cannot use system_to_component_index() ?

Best regards,

Sylvain Mathonnière

[Wolfgang Bangerth]

> Here *n_dofs = 320* and so *system_to_component_index()* does not like to have

> an argument i>=160 and throws an error.
> How come n_dofs=320 in this case ? I was expecting 2 interfaces * 2 nodes * 40
> components = 160 dofs, unless I am considering all the interfaces of the cell
> at onces in which case I have 2 interfaces * 4 nodes * 40 components = 320.
>

> *My main question is therefore:*

> How can I access the ith component of my shape function in face_worker like in
> the cell_worker if I cannot use system_to_component_index() ?

You are using FEInterfaceValues, which considers the shape functions that live
on the interface. This set of shape functions is the union of the shape
functions living on the two adjacent cells. Because you have a discontinuous
element, this set has size 2 * the number of dofs per cell, which gives you
exactly the 2*160=320 you observe.

Since you can't use FiniteElement::system_to_component_index on your index
'i', you need to find a different way. One approach is to compute a mapping
from this index 'i' to a pair '(here_or_there, index within cell)' and then
call system_to_component_index() on the second half of the pair.

Best
W.

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

[Sylvain Mathonnière]

Thank you for the answer. I have been looking at a way of doing what you suggested and I found the function FEInterfaceValues::interface_dof_to_dof_indices() which seems to be doing what I want.

I provide it with the interface index "i" and it returns a pair of local indices corresponding to the active and neighbouring cell. If then I query the first element  dof_index[0], I obtain the local index of the current cell, so I though of using it like this :

            //get the local indices of the current cell and neighbouring cell

            std::array<unsigned int, 2> dof_index = fe_iv.interface_dof_to_dof_indices(i);

 

           //use local index to obtain component of shape function.

           unsigned int component_i = fe_RTE.system_to_component_index(dof_index[0]).first;

Then with the debugger and printing some values I realised that dof_index[0] is kind of always below 160 as expected, but on a couple of occasion it goes to 22116; 89130 or 4294967295 and triggers the error of system_to_component_index(). Those values looks like global indices whereas interface_dof_to_dof_indices(i) should returns local_indices. I must be doing something horribly wrong but I could not find a tutorial that uses interface_dof_to_dof_indices(). Am I confused with how to use those functions ?

Best,

Sylvain

[Wolfgang Bangerth]

On 7/23/21 8:59 AM, Sylvain Mathonnière wrote:
> Thank you for the answer. I have been looking at a way of doing what you

> suggested and I found the function*FEInterfaceValues
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cbe067ec5c8ab4601e70308d94dea8304%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637626491881337201%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=hrnX9t7udM50eZvFVTiCBunh8ShGkow1iyI4AZq7udw%3D&reserved=0>::interface_dof_to_dof_indices()*

> which seems to be doing what I want.

Ah yes, I had forgotten about this function!

> I provide it with the interface index "i" and it returns a pair of local
> indices corresponding to the active and neighbouring cell. If then I query the

> first element *dof_index[0]*, I obtain the local index of the current cell, so

> I though of using it like this :
>

> *            //get the local indices of the current cell and neighbouring cell*
> **
> *            std::array<unsigned int, 2> dof_index =
> fe_iv.interface_dof_to_dof_indices(i);*
> **
> **
> **
> *           //use local index to obtain component of shape function.*
> **
> *           unsigned int component_i =
> fe_RTE.system_to_component_index(dof_index[0]).first;*

>
> Then with the debugger and printing some values I realised that dof_index[0]
> is kind of always below 160 as expected, but on a couple of occasion it goes
> to 22116; 89130 or 4294967295 and triggers the error of
> system_to_component_index(). Those values looks like global indices whereas
> interface_dof_to_dof_indices(i) should returns local_indices. I must be doing
> something horribly wrong but I could not find a tutorial that uses
> interface_dof_to_dof_indices(). Am I confused with how to use those functions ?

No, I think you're right that that is how the function is supposed to work.

4294967295 is numbers::invalid_dof_index, which is the case discussed in the
documentation where a DoF lives only one one side of the interface. The other
values (22116; 89130) shouldn't happen. If that's what you get, then there is
a bug somewhere.

Do you think you could construct a small testcase that illustrates what you have?

[Sylvain Mathonnière]

Oh, I think I know what those values appears and it is my bad ...

It is probably not 22116 but 22 and 116 and it is probably not 89130 but rather 89 and 130. I used something like std::cout << dof_index[0] << std::endl to get those values and sometimes they just get attached to each other like if the std::endl was skipped (not sure why, maybe because I have several threads ?). I did not realise it before reading your message but it actually makes sense because I was not getting an error for those values and my Assert was not triggered. I should have seen it, sorry.

As for the 4294967295 occurence, is it because I am dealing with an interface which is a boundary as well ? Isn't this case handled by the boundary_worker and not the face_worker (using step-12 notation) ?

Best,

Sylvain

[Corbin Foucart]

Sylvain,

Oh, I think I know what those values appears and it is my bad ...

It is probably not 22116 but 22 and 116 and it is probably not 89130 but rather 89 and 130. I used something like std::cout << dof_index[0] << std::endl to get those values and sometimes they just get attached to each other like if the std::endl was skipped (not sure why, maybe because I have several threads ?). I did not realise it before reading your message but it actually makes sense because I was not getting an error for those values and my Assert was not triggered. I should have seen it, sorry.

 

I ran into this exact issue a few weeks ago. For printing statements while debugging in assembly loops, check the signatures of functions like MeshWorker::run and MeshWorker::mesh_loop. You should be able to pass parameters that force non-multithreaded execution, and your print statements will make sense. I found this helpful for debugging.

Best,

Corbin

[Wolfgang Bangerth]

On 7/23/21 9:41 AM, Sylvain Mathonnière wrote:
>
> As for the 4294967295 occurence, is it because I am dealing with an interface

> which is a boundary as well ? Isn't this case handled by the *boundary_worker*
> and not the *face_worker* (using step-12 notation) ?

No, it's because you are using DG elements. I tried to clarify this here:
https://github.com/dealii/dealii/pull/12595

[Sylvain Mathonnière]

Thank you for the github link. In the end I fixed the issue by checking which value of the dof_index vector was not equal to numbers::invalid_unsigned_int. I added this line :

unsigned int component_i = (dof_index[0] != numbers::invalid_unsigned_int) ? 

                                    fe_RTE.system_to_component_index(dof_index[0]).first:

                                    fe_RTE.system_to_component_index(dof_index[1]).first;

Thank you Corbin as well. I will try that next time I am in need of it !

Best,

Sylvain M.

[Wolfgang Bangerth]

On 7/26/21 8:47 AM, Sylvain Mathonnière wrote:
> *
> *
> **
> *unsignedintcomponent_i = (dof_index[0] != numbers::invalid_unsigned_int) ? *
> *fe_RTE.system_to_component_index(dof_index[0]).first:*
> *fe_RTE.system_to_component_index(dof_index[1]).first;*

Yes, that's exactly the way to go.