Step-37 get_ function_values equivalent

[sheshu]

Hi,

I am trying to use matrix free methods based on step-37 for the elliptic equation in my system. The coefficient in this case depends on a variable "u" I get by solving an advection equation using DG methods. What is the best way to get these function values of u at the quadrature points for matrix free methods?

Thanks,
Sheshu

[Martin Kronbichler]

Dear Sheshu,

I am trying to use matrix free methods based on step-37 for the elliptic equation in my system. The coefficient in this case depends on a variable "u" I get by solving an advection equation using DG methods. What is the best way to get these function values of u at the quadrature points for matrix free methods?

What you need to do is to add the DoFHandler of the DG method to the reinit function of MatrixFree (check for the reinit functions with a std::vector<const DoFHandler<dim>*> argument and add one DoFHandler for the elliptic equation and one for the DG DoFHandler). Then, in an apply_local function, you simply create something like

  FEEvaluation<dim,fe_degree_dg,n_q_points_1d,n_components> dg_dofs(matrix_free, 1);

  FEEvaluation<dim,fe_degree,n_q_points_1d,n_components> elliptic(matrix_free, 0);

(assuming you put the elliptic DoFHandler into the zeroth slot and the DG DoFHandler in the first one). On a "macro cell", you call

  dg_dofs.reinit(cell);

  dg_dofs.read_dof_values(dg_solution_vector);

  dg_dofs.evaluate(true, false);

and on all quadrature points you can query dg_dofs.get_value(q) and do whatever you want with that (or similarly if you need the gradient...)

Best,

Martin

[sheshu]

Dear Martin,

Thanks a lot for your help. I tried following your reply but I get this error when I tried to use the reinit function from the MatrixFree class.
I am trying to call the reinit function defined at:
http://www.dealii.org/developer/doxygen/deal.II/classMatrixFree.html#a2e07eda8a73a7ed9f5e1f6ca0ec93959

As a simple case if I just add the following lines of code to step-48 and comment out everything after "make_grid_and_dofs" function in that program.

 std::vector<DoFHandler<dim>* > combined_dof_handler;
 combined_dof_handler.push_back(&dof_handler);
 combined_dof_handler.push_back(&dof_handler);

 std::vector<ConstraintMatrix* > combined_constraints;
 combined_constraints.push_back(&constraints);
 combined_constraints.push_back(&constraints);


    matrix_free_data.reinit (combined_dof_handler, combined_constraints,
                             quadrature, additional_data);
    //matrix_free_data.reinit (dof_handler, constraints,
    //                              quadrature, additional_data);

When I do this , I get the error
step-48.cc:500: error: no matching function for call to ‘dealii::MatrixFree<2, double>::reinit(std::vector<dealii::DoFHandler<2, 2>*, std::allocator<dealii::DoFHandler<2, 2>*> >&, std::vector<dealii::ConstraintMatrix*, std::allocator<dealii::ConstraintMatrix*> >&, dealii::QGaussLobatto<1>&, dealii::MatrixFree<2, double>::AdditionalData&)’

I am using deal.ii 7.2.0. Sorry if I am doing a silly mistake somewhere.

[Martin Kronbichler]

Dear Sheshu,

> std::vector<DoFHandler<dim>* > combined_dof_handler;
> combined_dof_handler.push_back(&dof_handler);
> combined_dof_handler.push_back(&dof_handler);

You need to use a vector of constant DoFHandlers - std::vector<const
DoFHandler<dim>* >. MatrixFree does not modify these data structures
internally and the constness is to signal it more explicitly. And
similarly for std::vector<const ConstraintMatrix*>. Can you try again?

Best,
Martin

[sheshu]

Dear Martin,

That works perfectly.

Thanks :)
Sheshu