FESystem with FE_values_view returns wrong dimensions

[Maxi Miller]

I created a FESystem containing three FE_Q-elements with <dim>-dimensions each, using the suggested approach in https://www.dealii.org/8.4.0/doxygen/deal.II/classFESystem.html:

   

 template <int dim>
    std::vector<const FiniteElement<dim>*> Diffusion<dim>::create_fe_list(const unsigned int polynomial_degree)
    {
        std::vector<const FiniteElement<dim>*> fe_list;
        fe_list.push_back(new FE_Q<dim>(polynomial_degree));
        return fe_list;
    }

    template <int dim>
    std::vector<unsigned int> Diffusion<dim>::create_fe_multiplicities()
    {
        std::vector<unsigned int> multiplicities;
        multiplicities.push_back(3);
        return multiplicities;
    }

    template <int dim>
    Diffusion<dim>::VectorElementDestroyer::VectorElementDestroyer(const std::vector<const FiniteElement<dim> *> &pointers) : data(pointers)
    {}

    template <int dim>
    Diffusion<dim>::VectorElementDestroyer::~VectorElementDestroyer()
    {
        for(size_t i = 0; i < data.size(); ++i)
            delete data[i];
    }template <int dim>
    std::vector<const FiniteElement<dim>*> Diffusion<dim>::create_fe_list(const unsigned int polynomial_degree)
    {
        std::vector<const FiniteElement<dim>*> fe_list;
        fe_list.push_back(new FE_Q<dim>(polynomial_degree));
        return fe_list;
    }

    template <int dim>
    std::vector<unsigned int> Diffusion<dim>::create_fe_multiplicities()
    {
        std::vector<unsigned int> multiplicities;
        multiplicities.push_back(3);
        return multiplicities;
    }

    template <int dim>
    Diffusion<dim>::VectorElementDestroyer::VectorElementDestroyer(const std::vector<const FiniteElement<dim> *> &pointers) : data(pointers)
    {}

    template <int dim>
    Diffusion<dim>::VectorElementDestroyer::~VectorElementDestroyer()
    {
        for(size_t i = 0; i < data.size(); ++i)
            delete data[i];
    }

    template <int dim>
    const std::vector<const FiniteElement<dim>*> &Diffusion<dim>::VectorElementDestroyer::get_data() const
    {
        return data;
    }

    finite_element(VectorElementDestroyer(create_fe_list(2)).get_data(), create_fe_multiplicities()),

    template <int dim>
    const std::vector<const FiniteElement<dim>*> &Diffusion<dim>::VectorElementDestroyer::get_data() const
    {
        return data;
    }

    finite_element(VectorElementDestroyer(create_fe_list(2)).get_data(), create_fe_multiplicities())

In addition I created in a later function corresponding views in order to access the different FE_Q-elements:

const FEValuesExtractors::Vector N_density(0);
const FEValuesExtractors::Vector E_temperature(dim);
const FEValuesExtractors::Vector L_temperature(2*dim);

and the corresponding code for DoFHandler and other parts:

system_matrix = 0.;
mass_matrix = 0.;

const QGauss<dim> quadrature_formula(fe_degree+1);

FEValues<dim> fe_values(finite_element, quadrature_formula,
                              update_values | update_gradients | update_JxW_values);


const unsigned int dofs_per_cell = finite_element.dofs_per_cell;
const unsigned int n_q_points    = quadrature_formula.size();

FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);

std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

std::vector<Tensor<2, dim>> gradient_N(dofs_per_cell), gradient_TE(dofs_per_cell), gradient_TL(dofs_per_cell);
std::vector<Tensor<1, dim>> value_N(dofs_per_cell), value_TE(dofs_per_cell), value_TL(dofs_per_cell);

for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
{
     cell_matrix = 0.;
          cell_mass_matrix = 0.;

          fe_values.reinit (cell);

          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
          {
              for(size_t i = 0; i < dofs_per_cell; ++i)
              {
                    gradient_N[i] = fe_values[N_density].gradient(i, q_point);
                    gradient_TE[i] = fe_values[E_temperature].gradient(i, q_point);
                    gradient_TL[i] = fe_values[L_temperature].gradient(i, q_point);
                    value_N[i] = fe_values[N_density].value(i, q_point);
                    value_TE[i] = fe_values[E_temperature].value(i, q_point);
                    value_TL[i] = fe_values[L_temperature].value(i, q_point);
              }

Compilation works fine for dim=3, but when executing I get the following error(s):

--------------------------------------------------------
An error occurred in line <3996> of file </opt/dealII/include/deal.II/fe/fe_values.h> in function
    const dealii::FEValuesViews::Vector<dim, spacedim>& dealii::FEValuesBase<dim, spacedim>::operator[](const dealii::FEValuesExtractors::Vector&) const [with int dim = 3; int spacedim = 3]
The violated condition was:
    vector.first_vector_component < fe_values_views_cache.vectors.size()
Additional information:
    Index 3 is not in the half-open range [0,1).

Stacktrace:
-----------
#0  ./ttm_test_2: Step52::Diffusion<3>::assemble_system(double)
#1  ./ttm_test_2: Step52::Diffusion<3>::run()
#2  ./ttm_test_2: main
--------------------------------------------------------

Where did I make my mistake?

[Wolfgang Bangerth]

>
> |
> template<intdim>
>
> std::vector<constFiniteElement<dim>*>Diffusion<dim>::create_fe_list(constunsignedintpolynomial_degree)
> {
> std::vector<constFiniteElement<dim>*>fe_list;
> fe_list.push_back(newFE_Q<dim>(polynomial_degree));
> returnfe_list;

So you have one scalar base element...


> template<intdim>
> std::vector<unsignedint>Diffusion<dim>::create_fe_multiplicities()

> {
> std::vector<unsignedint>multiplicities;
> multiplicities.push_back(3);

...with multiplicity three, for a total of three components. As a consequence...

> constFEValuesExtractors::VectorE_temperature(dim);
> constFEValuesExtractors::VectorL_temperature(2*dim);

...using these fill have to fail if dim>=2. Indeed, this happens:


> --------------------------------------------------------
> Anerror occurred inline <3996>of file
> </opt/dealII/include/deal.II/fe/fe_values.h>infunction
> constdealii::FEValuesViews::Vector<dim,spacedim>&dealii::FEValuesBase<dim,spacedim>::operator[](constdealii::FEValuesExtractors::Vector&)const[withintdim
> =3;intspacedim =3]
> Theviolated condition was:
> vector.first_vector_component <fe_values_views_cache.vectors.size()
> Additionalinformation:
> Index3isnotinthe half-open range [0,1).

Best
W.

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

[Maxi Miller]

Why do I have a scalar base element when initializing FE_Q<dim> with dim=3? I thought that I then have the finite element as 3d-element?

[Wolfgang Bangerth]

On 08/03/2017 08:47 AM, Maxi Miller wrote:
> Why do I have a scalar base element when initializing FE_Q<dim> with dim=3? I
> thought that I then have the finite element as 3d-element?

FE_Q is a scalar element, regardless of the dimension. It has only one vector
component.

[Maxi Miller]

I.e. I should rather choose FESystem<dim> a(FE_Q<dim>(2), 3*dim)?

[Wolfgang Bangerth]

On 08/03/2017 09:01 AM, 'Maxi Miller' via deal.II User Group wrote:
> I.e. I should rather choose FESystem<dim> a(FE_Q<dim>(2), 3*dim)?

If you want to describe three vector-fields, then that seems appropriate.