FESystem and FEInterfaceValues: how to get component

[Praveen C]

Hello

While assembling interface terms in a DG code with FESystem and FEInterfaceValues, I am not able to get the component which I would normally do with

for(i=0; i<fe.dofs_per_cell; ++i)
fe.system_to_component_index(i).first

Now since “i” comes from joint dof indices, above will not work.

Do you have a suggestion for this ?

I suppose I could use interface_dof_to_dof_indices()

https://dealii.org/current/doxygen/deal.II/classFEInterfaceValues.html#ad3e0203907da47ea2e312e94ec342b39

which gives the dof index on one cell and then do fe.system_to_component_index(i).first.

But this returns a tuple, and I have to find one of them which is valid, which seems ugly.

Thanks
praveen

[Praveen C]

I am doing the following in my face worker, which seems to work.

best

praveen

   const unsigned int n_cell_dofs = fe_face_values.get_fe().n_dofs_per_cell();

   const unsigned int n_face_dofs = fe_face_values.n_current_interface_dofs();

   for(unsigned int q=0; q<n_q_points; ++q)

   {

      Vector<double> num_flux(nvar);

      PDE::numerical_flux(param->flux_type, 

                          left_state[q], 

                          right_state[q], 

                          q_points[q], 

                          fe_face_values.normal(q),

                          num_flux);

      for (unsigned int i = 0; i < n_face_dofs; ++i)

      {

         unsigned int ii = (i < n_cell_dofs) ? i : i - n_cell_dofs;

         const auto c = fe_face_values.get_fe().system_to_component_index(ii).first;

         cell_rhs(i) -= num_flux[c] *

                        fe_face_values.jump_in_shape_values(i, q, c) *

                        fe_face_values.JxW(q);

[Wolfgang Bangerth]

On 9/14/24 08:04, Praveen C wrote:
>          unsigned int ii = (i < n_cell_dofs) ? i : i - n_cell_dofs;
>          const auto c =
> fe_face_values.get_fe().system_to_component_index(ii).first;

Praveen:
this makes the assumption that you first have the DoFs from one side, and then
the ones from the other side. This may be true in some cases (including,
perhaps, yours) but it is not one that you should rely on unless you have a
way of asserting that your assumption is correct.

I don't think we have a good way of querying this information at the moment.
You might have to ask for the DoF indices of both cells, and then of the
FEInterfaceValues object, and compare these indices so that when you need
information for the FEInterfaceValues object, you know which of the two cells
to ask and how.

I think a better way would actually be to just implement this sort of
functionality in FEInterfaceValues to begin with. The class is not very
difficult to understand, and it shouldn't be very difficult to implement what
you need. The usual applies: We'd love to have a patch that we can incorporate
into the library!

Best
W.

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

[Praveen C]

The simplest solution seems to be to add a function like this.

// This is available only after call to reinit.

unsigned int

FEInterfaceValues::get_component(const unsigned int interface_dof_index) const

{

   const auto dof_pair = dofmap[interface_dof_index];

   if(dof_pair[0] != numbers::invalid_unsigned_int)

   {

      return get_fe().system_to_component_index(dof_pair[0]).first;

   }

   else

   {

      return get_fe().system_to_component_index(dof_pair[1]).first;

   }

}

Or should the component info be built during reinit and stored in the class ?

best

praveen

--
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 the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/80531e2b-863d-45b5-a0fa-68859f69181e%40colostate.edu.

[Praveen C]

It should look like this so it works in hp case.

// This is available only after call to reinit.

unsigned int

FEInterfaceValues::get_component(const unsigned int interface_dof_index) const

{

   const auto dof_pair = dofmap[interface_dof_index];

   if(dof_pair[0] != numbers::invalid_unsigned_int)

   {

      return fe_face_values.get_fe().system_to_component_index(dof_pair[0]).first;

   }

   else

   {

      return fe_face_values_neighbor.get_fe().system_to_component_index(dof_pair[1]).first;

   }

}

best

praveen

[Wolfgang Bangerth]

On 9/15/24 23:07, Praveen C wrote:
> **

>
> The simplest solution seems to be to add a function like this.
>
> // This is available only after call to reinit.
> unsigned int
> FEInterfaceValues::get_component(const unsigned int interface_dof_index) const
> {
>    const auto dof_pair = dofmap[interface_dof_index];
>    if(dof_pair[0] != numbers::invalid_unsigned_int)
>    {
>       return get_fe().system_to_component_index(dof_pair[0]).first;
>    }
>    else
>    {
>       return get_fe().system_to_component_index(dof_pair[1]).first;
>    }
> }
>
> Or should the component info be built during reinit and stored in the class ?

This (or your later variation) works just fine. I don't think you'd save
anything by computing it in reinit().

I would add an assertion at the top that interface_dof_index is within the
allowed range.