Error calculating shape gradient

[A.Z Ihsan]

Hi all, 

i am getting this error while calculating the local cell matrix for my hp fem application

dealii::Tensor<1, spacedim> dealii::FEValuesBase<dim, spacedim>::shape_grad_component(unsigned int, unsigned int, unsigned int) const [with int dim = 3; int spacedim = 3]

The violated condition was: 

    component < fe->n_components()

Additional information: 

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

Stacktrace:

-----------

#0  ../build/local: dealii::FEValuesBase<3, 3>::shape_grad_component(unsigned int, unsigned int, unsigned int) const

#1  ../build/local: dealii::SymmetricTensor<2, 3, double> LinearElasticity::get_strain<3>(dealii::FEValues<3, 3> const&, unsigned int, unsigned int)

#2  ../build/local: LinearElasticity::HPSolver<3, dealii::Vector<double> >::assemble_cell_matrix(dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::hp::DoFHandler<3, 3>, false> > const&, dealii::FullMatrix<double>&, dealii::hp::FEValues<3, 3>&)

#3  ../build/local: LinearElasticity::HPSerialSolver<3>::assemble_linear_system(LinearElasticity::HPSerialSolver<3>::LinearSystem&)

#4  ../build/local: LinearElasticity::HPSerialSolver<3>::solve_problem()

 

The idea is similar to the tutorial step-27, but here i use the symmetric tensor for returning strain. Can someone explain to me what the error is and how to resolve it? 

BR, 

Ihsan

[Simon Sticko]

Hi,
from the error message you can see that the element you are using only has 1 component. You get an error because you are trying to access component 1, which doesn't exist. Since your element should have dim components, there is likely something wrong with how you create your element. It should probably be similar as in step-8.

Best,
Simon

[A.Z Ihsan]

Hi Simon, 

i followed the tutorial step-27 about hp-fem, using the similar constructor with dim = 3

for (unsigned int deg=1; deg <= max_fe_degree; ++deg)

      {

        fe_collection.push_back(FE_Q<dim>(deg));

        quadrature_collection.push_back(QGauss<dim>(deg+1));

        face_quadrature_collection.push_back(QGauss<dim-1>(deg+1));

      }

I tried to look into the fe_collection.n_components(), it yields 1.

is that what you meant?

BR, 

Ihsan

[Simon Sticko]

Hi,
Yes. Here, you add a number of scalar elements of different degree:

> for (unsigned int deg=1; deg <= max_fe_degree; ++deg)
> {
> fe_collection.push_back(FE_Q<dim>(deg));
> quadrature_collection.push_back(QGauss<dim>(deg+1));
> face_quadrature_collection.push_back(QGauss<dim-1>(deg+1));

> }

which makes sense for a scalar problem, like the Poisson equation.

But if you solve for displacement in a material the solution has dim components:

u = (u_x, u_y, u_z)

so you want to use a FESystem element with dim components:

fe_collection.push_back(FESystem<dim>(FE_Q<dim>(deg), dim));

Best,
Simon

[A.Z Ihsan]

Great, it works now. 

THanks Simon. 

BR, 

Ihsan