Matrix-Free and FE_Nothing

[Thomas Cheetham]

Hi all,

I am trying to implement a MatrixFree code which uses FE_Q and FE_Nothing elements. Evaluating the matrix-free operator works correctly without the FE_Nothing elements but when I add them the code the operator tries to access indices which are out of range in the read_dof_values function. I was wondering if anyone has any tips for applying this correctly?

Here is my current local_apply function:

template <int dim, int degree_v, typename SignedDistanceType, typename number>

void MassOperator<dim, degree_v, SignedDistanceType, number>::local_apply(

    const dealii::MatrixFree<dim, number>                 &data,

    LinearAlgebra::distributed::BlockVector<number>       &dst,

    const LinearAlgebra::distributed::BlockVector<number> &src,

    const std::pair<unsigned int, unsigned int>           &cell_range) const

{

    using vector_t = VectorizedArray<number>;

    FEEvaluation<dim, degree_v, degree_v*2, dim, number>   velocity(data, 0);

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

    {

        //AssertDimension(phase_field_evaluations.size(0), data.n_cell_batches());

        //AssertDimension(phase_field_evaluations.size(1), velocity.n_q_points);

        velocity.reinit(cell);

        //std::cout << velocity.get_active_fe_index() << std::endl;

        if (velocity.get_active_fe_index()==0){

            // Dof are active and not FENothing

            velocity.read_dof_values(src.block(0));

   

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

            {

                //velocity.submit_value(phase_field_evaluations(cell, q)*velocity.get_value(q),q);

                velocity.submit_value(velocity.get_value(q),q);

            }

            velocity.integrate(EvaluationFlags::values);

            velocity.distribute_local_to_global(dst.block(0));

           

        }

    }

}

Regards,

Thomas

[Bruno Turcksin]

Hello Thomas,

If you use FE_Nothing, you don't want to loop over the cell_range. Instead you create a subrange of cells using create_cell_subrange_hp_by_index() and then loop over the cells in the subrange.

Best,

Bruno

[Thomas Cheetham]

Hi Bruno,

Thanks for your quick reply, that worked perfectly.

Many thanks,

Thomas