Is MatrixFreeOperators::LaplaceOperator::local_diagonal_cell() something wrong?

[llf m]

Dear all,

The function MatrixFreeOperators::LaplaceOperator::local_diagonal_cell() in matrixfree/operators.h line 2135-2153

  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(

    data, this->selected_rows[0]);

  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)

    {

      phi.reinit(cell);

      VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];

      for (unsigned int i = 0; i < phi.dofs_per_component; ++i)

        {

          for (unsigned int j = 0; j < phi.dofs_per_component; ++j)

            phi.begin_dof_values()[j] = VectorizedArrayType();

          phi.begin_dof_values()[i] = 1.;

          do_operation_on_cell(phi, cell);

          local_diagonal_vector[i] = phi.begin_dof_values()[i];

        }

      for (unsigned int i = 0; i < phi.dofs_per_component; ++i)

        for (unsigned int c = 0; c < phi.n_components; ++c)

          phi.begin_dof_values()[i + c * phi.dofs_per_component] =

            local_diagonal_vector[i];

      phi.distribute_local_to_global(dst);

In the first loop inside the cell loop index i go through from 0 to phi.dofs_per_component (highlighted) as well as the j index,

but the function phi.begin_dof_values()return a pointer whose index go throught from 0 the phi.dofs_per_cell, that is said

in n_component > 1 case the index after phi.dofs_per_component never been touched( the total is dofs_per_compoenent*n_component).

I have found an old discussion in https://groups.google.com/forum/#!searchin/dealii/MatrixFreeOperators$3A$3ALaplaceOperator%7Csort:date/dealii/k133wjcKM70/Jw-reYrsAwAJ\

but things change now.

Am I misunderstand here?

best,

m.

[Martin Kronbichler]

In the first loop inside the cell loop index i go through from 0 to phi.dofs_per_component (highlighted) as well as the j index,

but the function phi.begin_dof_values()return a pointer whose index go throught from 0 the phi.dofs_per_cell, that is said

in n_component > 1 case the index after phi.dofs_per_component never been touched( the total is dofs_per_compoenent*n_component).

I have found an old discussion in https://groups.google.com/forum/#!searchin/dealii/MatrixFreeOperators$3A$3ALaplaceOperator%7Csort:date/dealii/k133wjcKM70/Jw-reYrsAwAJ\

but things change now.

The code here should be correct, because we are looking at a constant-coefficient Laplacian here, where all entries to the diagonal are equal for the various components. So even though the local computation does indeed compute all contributions (but with some spurious dummy entries in the components higher than zero; which eventually are never used), but only get the zero-th component of it and duplicate it over all entries. This is done to save some work, given that we know that the entry is the same; in case you do an operator that couples between the components, one of course has to go over all phi.dofs_per_cell components both when zeroing entries and in the loop over the columns i.

Best,
Martin

[Martin Kronbichler]

I should add that you can set a variable coefficient, but it is still scalar and thus the same to all components, so the code is correct. It would be different if we allowed tensor-valued coefficients (which we do not at the moment); then, one would really need to go through all dofs_per_cell entries.

[llf m]

Dear Martin,

I get the meaning, you are right. 

And in this case 

if

VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];

change to

VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_component];

should also work, right?
在 2019年11月15日星期五 UTC+8下午3:27:44，Martin Kronbichler写道：