Matrix-free: Continuous elements and boundary (face) integrals

[Michał Wichrowski]

Deal all,

I want to weakly impose boundary conditions on, let's say, Laplace problem (in fact I'm working on Stokes).  I'm considering [FE_Q<2>(2)^2]   ( that is   FESystem<2>[FE_Q<2>(2)^2] + FE_Q<2>(1) for Stokes).

As a MatrixFree loop I'm using:

  this->data->loop(&SymGradMass::local_apply, &SymGradMass::apply_face,

                   &SymGradMass::apply_boundary,

   this, dst, src,

                   /*zero_dst =*/false,

                   MatrixFree<dim, Number>::DataAccessOnFaces::gradients,

                   MatrixFree<dim, Number>::DataAccessOnFaces::gradients);

where apply_face function is empty (only local_apply and apply_boundary perform integration ).

 After initializing MatrixFree object with following data:

    typename MatrixFree<dim, Number>::AdditionalData additional_data;

    additional_data.tasks_parallel_scheme =

    MatrixFree<dim,Number>::AdditionalData::none;

    additional_data.mapping_update_flags =

      (update_gradients | update_JxW_values | update_quadrature_points);

    additional_data.mapping_update_flags_boundary_faces =

        (update_gradients | update_JxW_values | update_normal_vectors |

         update_quadrature_points);

    additional_data.mapping_update_flags_inner_faces =

        (update_default);

I get end error (at the end of the email).  It looks like the MatrixFree is trying to compute coupling of DoF over face but on the other side that information is not provided (by FE?). I suppose this only requires simple fix:  in case if 

mapping_update_flags_inner_faces =    update_default

proceed as in standard case (no boundary integrals) - the boundary integrals does not generate additional DoF coupling.

Best,

Michał Wichrowski

An error occurred in line <786> of file </home/mwichro/lib/dealii/include/deal.II/base/partitioner.h> in function

    dealii::types::global_dof_index dealii::Utilities::MPI::Partitioner::local_to_global(unsigned int) const

The violated condition was: 

    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(local_index), decltype(local_size() + n_ghost_indices_data)>::type)>:: type>(local_index) < static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(local_index), decltype(local_size() + n_ghost_indices_data)>::type)>:: type>(local_size() + n_ghost_indices_data)

Additional information: 

    Index 4294967295 is not in the half-open range [0,12832).

Stacktrace:

-----------

#0  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: dealii::Utilities::MPI::Partitioner::local_to_global(unsigned int) const

#1  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: void dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 4> >::initialize_indices<double>(std::vector<dealii::AffineConstraints<double> const*, std::allocator<dealii::AffineConstraints<double> const*> > const&, std::vector<dealii::IndexSet, std::allocator<dealii::IndexSet> > const&, dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 4> >::AdditionalData const&)

#2  /home/mwichro/lib/deal.II/lib/libdeal_II.g.so.9.2.0-pre: void dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 4> >::internal_reinit<double>(dealii::Mapping<2, 2> const&, std::vector<dealii::DoFHandler<2, 2> const*, std::allocator<dealii::DoFHandler<2, 2> const*> > const&, std::vector<dealii::AffineConstraints<double> const*, std::allocator<dealii::AffineConstraints<double> const*> > const&, std::vector<dealii::IndexSet, std::allocator<dealii::IndexSet> > const&, std::vector<dealii::hp::QCollection<1>, std::allocator<dealii::hp::QCollection<1> > > const&, dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 4> >::AdditionalData const&)

#3  ./fsi_matrix_free_2d.g: void dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 4> >::reinit<dealii::DoFHandler<2, 2>, dealii::QGauss<1>, double>(dealii::Mapping<2, 2> const&, std::vector<dealii::DoFHandler<2, 2> const*, std::allocator<dealii::DoFHandler<2, 2> const*> > const&, std::vector<dealii::AffineConstraints<double> const*, std::allocator<dealii::AffineConstraints<double> const*> > const&, dealii::QGauss<1> const&, dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 4> >::AdditionalData const&)

#4  ./fsi_matrix_free_2d.g: StokesSolver<2, 2>::setup_system()

#5  ./fsi_matrix_free_2d.g: StokesSolver<2, 2>::initialize(std::shared_ptr<dealii::Mapping<2, 2> >&)

#6  ./fsi_matrix_free_2d.g: FSIInterface<2, 2>::FSIInterface(RunTimeParameters const&)

#7  ./fsi_matrix_free_2d.g: main

--------------------------------------------------------

[Michał Wichrowski]

I tried to fix that by myself, but I could not locate the source code of MatrixFree::internal_reinit

[Bruno Turcksin]

On Wednesday, March 11, 2020 at 7:36:31 AM UTC-4, Michał Wichrowski wrote:

I tried to fix that by myself, but I could not locate the source code of MatrixFree::internal_reinit

Here it is: https://github.com/dealii/dealii/blob/master/include/deal.II/matrix_free/matrix_free.templates.h#L352-L681 Note that they are two different internal_reinit functions.

Best,

Bruno