Ordering using multi-components

[Joss G.]

Dear all, 

It is not very clear to me what is the order followed when using multi-components. As an example to put my question more clear I cann point to step-29 where there are 2 components used.

My common sense would be that the system_matrix global operator would have first real and then the imaginar components. However, printing the matrix and looking at it, I can see the values are somehow intercalated (rows and columns).

Could someone please clarify this?

I would also like to ask how can I get the mass matrix of one of the components. Using the following for each cell and locating for the global matrix I get the two components (again not sure how the order is...)

Mass_matrix(i,j)+=((fe_values.shape_value(i, q_point) * fe_values.shape_value(j, q_point)) * fe_values.JxW(q_point));

On the other hand, using the following, does not allow me to use the local_dof_indices function, as using only one component, the dof is half:

Mass_matrix(i,j)+=((fe_values.shape_value_component(i, q_point,0) * fe_values.shape_value(j, q_point,0)) * fe_values.JxW(q_point));

What is the best way to get this?

Thank you very much for your time.

[Wolfgang Bangerth]

Joss:

> It is not very clear to me what is the order followed when using
> multi-components. As an example to put my question more clear I cann
> point to step-29 where there are 2 components used.
>
> My common sense would be that the system_matrix global operator would
> have first real and then the imaginar components. However, printing the
> matrix and looking at it, I can see the values are somehow intercalated
> (rows and columns).

> *Could someone please clarify this?*

The way we typically like to think about this is that the DoFHandler
chooses "some" order that simply is what it is. If you want something
specific, you need to say so explicitly. In your case, if you want all
of the real components first and then all of the imaginary components,
then you need to say that that is what you want, and the way to do that
is DoFRenumbering::component_wise(). This function is used in any number
of tutorial programs, see for example step-20 and step-22.

> I would also like to ask how can I get the mass matrix of one of the
> components. Using the following for each cell and locating for the
> global matrix I get the two components (again not sure how the order is...)

> /Mass_matrix(i,j)+=((fe_values.shape_value(i, q_point) *
> fe_values.shape_value(j, q_point)) * fe_values.JxW(q_point));/

This results in a different matrix than you want: If the matrix was
structured by component (i.e., it is a 2x2 block matrix), then each of
the four blocks will end up with a mass matrix. That's likely not what
you intended.

On the other hand...

>
> On the other hand, using the following, does not allow me to use the
> local_dof_indices function, as using only one component, the dof is half:

> /Mass_matrix(i,j)+=((fe_values.shape_value_component(i, q_point,0) *
> fe_values.shape_value(j, q_point,0)) * fe_values.JxW(q_point));/

...this will only put a mass matrix into the top left block.

The general strategy to working with all multi-component problems is to
use the "extractors". See again step-20, step-22, and the documentation
module about vector-valued problems.

Best
W.


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

[Joss G.]

Dear Wolfgang Bangerth,

Thank you very much for your answer and for pointing to the relevant steps examples.

I would like to be sure that I am doing the right thing. First, I used DoFRenumbering::component_wise(dof_handler) and got the desired order, thank you. Then I use the extractors

const FEValuesExtractors::Scalar real_component(0);

const FEValuesExtractors::Scalar imag_component(1);

Is it correct the following approach to get the global mass matrix of one component? I get a big 2x2 block matrix with non-zero values on the left top block. Is this corresponding to the global mass matrix of one of the components?

Mass_matrix_tmp(i,j)+=((fe_values[real_component].value(i, q_point) * fe_values[real_component].value(j, q_point)) * fe_values.JxW(q_point));

Thank you again.

Best regards,

[Wolfgang Bangerth]

On 11/9/21 8:24 AM, Joss G. wrote:
> *Is it correct the following approach to get the global mass matrix of
> one component? *I get a big 2x2 block matrix with non-zero values on the

> left top block. Is this corresponding to the global mass matrix of one
> of the components?

> /Mass_matrix_tmp(i,j)+=((fe_values[real_component].value(i, q_point) *
> fe_values[real_component].value(j, q_point)) * fe_values.JxW(q_point));/

Yes, this will work.

[Joss G.]

Dear Wolfgang Bangerth,

Thank you for your clarification.

it is still not very clear to me why I am not getting the same mass matrix when doing this two different approaches when, as far as I understand should be the same. In order to get a 2x2 block matrix containing the mass matrix for the top left and bottom right I do the following:

Mass_matrix_tmp(i,j)+=((fe_values[real_component].value(i, q_point) * fe_values[real_component].value(j, q_point)) * fe_values.JxW(q_point));

Mass_matrix_tmp(i,j)+=((fe_values[imag_component].value(i, q_point) * fe_values[imag_component].value(j, q_point)) * fe_values.JxW(q_point));

which gives different results compared to:

if (fe.system_to_component_index(i).first == fe.system_to_component_index(j).first)

Mass_matrix_tmp(i,j)+= ((fe_values.shape_value(i, q_point) * fe_values.shape_value(j, q_point)) * fe_values.JxW(q_point));

Could you please clarify why this two approaches do not give the same results?

Thank you again

Regards,

[Wolfgang Bangerth]

> it is still not very clear to me why I am not getting the same mass
> matrix when doing this two different approaches when, as far as I
> understand should be the same. In order to get a 2x2 block matrix
> containing the mass matrix for the top left and bottom right I do the
> following:

> /Mass_matrix_tmp(i,j)+=((fe_values[real_component].value(i, q_point) *
> fe_values[real_component].value(j, q_point)) * fe_values.JxW(q_point));/

> /Mass_matrix_tmp(i,j)+=((fe_values[imag_component].value(i, q_point) *
> fe_values[imag_component].value(j, q_point)) * fe_values.JxW(q_point));/

>
> which gives different results compared to:

> /if (fe.system_to_component_index(i).first ==
> fe.system_to_component_index(j).first)/
> /Mass_matrix_tmp(i,j)+= ((fe_values.shape_value(i, q_point) *
> fe_values.shape_value(j, q_point)) * fe_values.JxW(q_point));/
> /
> /

> Could you please clarify why this two approaches do not give the same
> results?

It looks to me like these two methods should give the same matrix. In
which way are they different?