Extracting Shape Laplacian

[ky...@math.uh.edu]

Hello all,

I am working on implementing some stabilization methods for the Navier-Stokes equations, these techniques often require the use of the laplacian of the shape functions. I know I get the hessian of the k'th shape function at quadratue point q via the following

Tensor<3, dim> hessian_phi = fe_values[velocity_extractor].hessian(k, q);

My question is understanding what does this return, and how to extract the laplacian from it. Does hessian_phi[ i ][ j ][ k ] = \frac{  \partial \varphi_{ i }   } {  \partial x_{ j } \partial x_{ k }   }?

If this is the case, is there a way to contract over the last two components to result in a rank 1 tensor that is the laplacian? My idea was the following

std::vector< Tensor<1, dim> > laplace_phi_u (dofs_per_cell);
Tensor<3, dim> hessian_phi = fe_values[velocity_extractor].hessian(k, q);
for( int i=0; i<dim; ++i )
 laplace_phi_u[ k ][ i ] = trace( hessian_phi[ i ] )

This doesn't work, because at runtime it tells me that only 0 assignment is allowed.


I've also seen that I can use shape_hessian_component since I am using standard Q2-Q1 elements,

int component_i = fe.system_to_component( k ).first;
Tensor<2, dim> hessian_phi = fe_values.shape_hessian_component( k, q, component_i );

In this case I can get the laplacian of the i'th component of the k'th shape function at quadrature point q by trace( hessian_phi ). But this will result in a rank 0 tensor. It would be much more convenient if I could convert this back into a rank 1 tensor as the laplacian of the k'th shape function. Again though I run into the same problem as above with only 0 assignment is allowed.

Thanks

[ky...@math.uh.edu]

I found that my idea above was not what I had typed in my code, I am no longer receiving the run time error. However I would still like some verification that the indices are as I think they are.

[Wolfgang Bangerth]

On 08/09/2017 03:46 PM, ky...@math.uh.edu wrote:
>
> Hello all,
>
> I am working on implementing some stabilization methods for the Navier-Stokes
> equations, these techniques often require the use of the laplacian of the
> shape functions. I know I get the hessian of the k'th shape function at
> quadratue point q via the following
>
> |
> Tensor<3,dim>hessian_phi =fe_values[velocity_extractor].hessian(k,q);
> |
>
> My question is understanding what does this return, and how to extract the
> laplacian from it. Does hessian_phi[ i ][ j ][ k ] = \frac{ \partial
> \varphi_{ i } } { \partial x_{ j } \partial x_{ k } }?

I don't recall, but it should be documented somewhere with
FEValuesViews::Vector IIRC.

> If this is the case, is there a way to contract over the last two components
> to result in a rank 1 tensor that is the laplacian?

Not as you are trying, but you can easily write the summation over the last
two indices by hand, of course.

> I've also seen that I can use shape_hessian_component since I am using
> standard Q2-Q1 elements,
>
> |
> intcomponent_i =fe.system_to_component(k ).first;

> Tensor<2,dim>hessian_phi =fe_values.shape_hessian_component(k,q,component_i );

> |
>
> In this case I can get the laplacian of the i'th component of the k'th shape
> function at quadrature point q by trace( hessian_phi ). But this will result
> in a rank 0 tensor.

Yes, it is the trace of the i'th component. You can then construct a
Tensor<1,dim> traces;
for (i=0...dim)
{
traces_phi = fe_values.shape_hessian_component (...);
traces[i] = trace(hessian_phi)
}

Best
W.

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

[ky...@math.uh.edu]

Thanks,

I did try looking at the documentation for FEValuesViews::Vector< dim, spacedim >::hessian(const unsigned int shape_function, const unsigned int q_point)  before coming here, and all it provides is

"Return the Hessian (the tensor of rank 2 of all second derivatives) of the vector components selected by this view, for the shape function and quadrature point selected by the arguments."

I looked in fe_values.h to see if it was commented in the implementation and all that I see there is

  hessian_type return_value;
  for (unsigned int d=0; d<dim; ++d)
    if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
       return_value[d]
         = fe_values->finite_element_output.shape_hessians[shape_function_data[shape_function].row_index[d]][q_point];

which I'm not quite sure how to parse. The comments in the header file only say "same as for the scalar case except that we have one more index." (the added index is d).

This is why I brought my question here.

[Wolfgang Bangerth]

On 08/09/2017 05:40 PM, ky...@math.uh.edu wrote:
>
> I did try looking at the documentation for FEValuesViews::Vector

> <http://www.dealii.org/8.5.0/doxygen/deal.II/classFEValuesViews_1_1Vector.html><

> dim, spacedim >::hessian(const unsigned int shape_function, const unsigned int
> q_point) before coming here, and all it provides is
>
> "Return the Hessian (the tensor of rank 2 of all second derivatives) of the
> vector components selected by this view, for the shape function and quadrature
> point selected by the arguments."

Yes, the FEValuesViews::Vector::hessian(i,q) function returns the Hessian of
the i'th shape function at the q'th quadrature point as a matrix. I think that
ought to be clear.

The get_function_hessians() function returns a vector of tensors, which are
the hessians of the solution at the quadrature points.