help with get_function_divergences for a vector solution

[Jane Lee]

Dear all, 

I am hoping someone can shed a bit more light on a problem I am having. The background of the system of equations i am solving: I have 4 DAEs (for pf, pr, vr, vf) and 2 time dependent equations (phi, t) as a system of 6 equations. I am solving the 4daes first then using the solutions into the time dependent equations (which need the solutions in the DAEs), and I am using these two solutions to rerun the 4 DAEs and so on. 

The right hand side of the weak forms of the DAEs all require information about each other, so at the initial timestep I am providing a guess, but from then on, i am using get_function_values or get_function_gradients/divergences. I have been using these fine for most of my variables, but i am having an issue with one of them. 

My 4 DAEs look basically like the Stokes' system, with 1 DAEs in mixed formulation (very similar to Darcy's) and another non-differential algebraic equation for vf. Now, for this equation/variable, I am using:

        FESystem<dim> vf_fe;

        DoFHandler<dim>     vf_dof_handler;

        ConstraintMatrix    vf_constraints;

        SparsityPattern      vf_sparsity_pattern;

        SparseMatrix<double> vf_system_matrix;

        SparseMatrix<double> vf_mass_matrix;

        Vector<double>       vf_solution;

        Vector<double>       vf_initial_solution;

        Vector<double>       vf_system_rhs;

and  vf_fe (FE_Q<dim>(pf_degree+1), dim) as vf is vector valued. 

For the others, like the example codes (step-20/22/31), i am using BlockVectors etc. eg,

        rock_fe (FE_Q<dim>(pr_degree+1), dim,

    FE_Q<dim>(pr_degree), 1)  

with BlockSparsityPattern      rock_sparsity_pattern;

        BlockSparseMatrix<double> rock_system_matrix;

        BlockVector<double> rock_solution;

        BlockVector<double> rock_system_rhs; 

and so on. 

I am of course using const FEValuesExtractors::Vector velocities (0); etc for most of the things I need. now out of the following, in the loop over all the cells: I have

                 pf_fe_values[pressure].get_function_values (pf_solution, pf_values);

                pf_fe_values[pressure].get_function_gradients (pf_solution, grad_pf_values);

                phi_fe_values.get_function_values (phi_solution, phi_values);

                phi_fe_values.get_function_gradients (phi_solution, grad_phi_values);

                vf_fe_values.get_function_values (vf_solution, vf_values);

                vf_fe_values.get_function_divergences (vf_solution, div_vf_values);

and the LAST one for get_function_divergences does not seem to work. It is telling me there is no member function called get_function_divergeneces. I am a little confused as elsewhere i use vr_fe_values[velocities].get_function_divergences (rock_solution, div_vr_values); and this works fine. 

I have also tried vf_fe_values[velocities].get_function_divergences (vf_solution, div_vf_values); - now this does not give me the same error, the program runs but i get:

An error occurred in line <2364> of file </scratch/leej/dealnew/dealii-8.5.0/source/fe/fe_values.cc> in function

    void dealii::internal::do_function_values(const Number*, const dealii::Table<2, double>&, const dealii::FiniteElement<dim, spacedim>&, const std::vector<unsigned int>&, dealii::VectorSlice<std::vector<VectorType> >&, bool, unsigned int) [with int dim = 2; int spacedim = 2; VectorType = dealii::Vector<double>; Number = double]

The violated condition was: 

    (values[i].size()) == (result_components)

Additional information: 

    Dimension 0 not equal to 2.

I also get the same no member function call when used with the get_function_values on the same variable with FEValuesExtractors and adding [velocities]

I would be grateful if someone could shed a bit more light on why I am getting these errors. I particularly do not understand why it is giving me a dimension error either. It seems like quite a trivial thing but it seems not? I wondered whether it is because it is not aware that the solution is a vector, but I'm pretty sure the fact that I have an FESystem sorts this out?

As further information, note that i am using

std::vector<Vector<double> >        vf_values (n_q_points);

std::vector<double>   div_vf_values (n_q_points);

I have also looped over the right cells etc using

typename DoFHandler<dim>::active_cell_iterator

            vf_cell = vf_dof_handler.begin_active();

as well as intialising correctly. 

Before this assmebly bit is run, the initial vf_solution is just interpolated with 

VectorTools::interpolate(vf_dof_handler, vfInitialFunction<dim>(), vf_solution);

where it is just zero

    template <int dim>

    double

    vfInitialFunction<dim>::value (const Point<dim>  &p,

                                   const unsigned int component) const

    {

        if (component == 0)

            return 0.0;

        else if (component == 1)

            return 0.0*p[0];

        return 0.0;

    }

    

Many thanks for your help in advance. 

[Wolfgang Bangerth]

> I am of course using const FEValuesExtractors::Vector velocities (0); etc for
> most of the things I need. now out of the following, in the loop over all the
> cells: I have
> pf_fe_values[pressure].get_function_values (pf_solution, pf_values);
> pf_fe_values[pressure].get_function_gradients (pf_solution,
> grad_pf_values);
> phi_fe_values.get_function_values (phi_solution, phi_values);
> phi_fe_values.get_function_gradients (phi_solution, grad_phi_values);
> vf_fe_values.get_function_values (vf_solution, vf_values);
> vf_fe_values.get_function_divergences (vf_solution, div_vf_values);
> and the LAST one for get_function_divergences does not seem to work. It is
> telling me there is no member function called get_function_divergeneces. I
> am

This is correct -- FEValues has functions such as get_function_values and
get_function_gradients, but not get_function_divergences. That's because
FEValues corresponds to a general collection of vector components. But the
divergence is an operation that is only defined if there are exactly dim
vector components, which FEValues can not assume to be the case.

But if you use a FEValuesExtractors argument and say
vr_fe_values[velocities]
then you are selecting a subset of exactly dim components, and for those the
divergence operation is well defined. Consequently...

> a little confused as elsewhere i
> use vr_fe_values[velocities].get_function_divergences (rock_solution,
> div_vr_values); and this works fine.

...this works.

> I have also tried vf_fe_values[velocities].get_function_divergences
> (vf_solution, div_vf_values); - now this does not give me the same error, the
> program runs but i get:
> An error occurred in line <2364> of file
> </scratch/leej/dealnew/dealii-8.5.0/source/fe/fe_values.cc> in function
> void dealii::internal::do_function_values(const Number*, const
> dealii::Table<2, double>&, const dealii::FiniteElement<dim, spacedim>&, const
> std::vector<unsigned int>&, dealii::VectorSlice<std::vector<VectorType> >&,
> bool, unsigned int) [with int dim = 2; int spacedim = 2; VectorType =
> dealii::Vector<double>; Number = double]
> The violated condition was:
> (values[i].size()) == (result_components)
> Additional information:
> Dimension 0 not equal to 2.

This is a completely separate issue. I don't know how exactly you got there,
but it seems to me that the output vector you passed to the function has the
wrong size. Did you size it correctly?

> I also get the same no member function call when used with the
> get_function_values on the same variable with FEValuesExtractors and adding
> [velocities]

That is again a separate issue -- the VectorView::get_function_values function
just takes different argument types.

Best
W.

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