On 8/13/19 5:34 PM, Shiva Rudraraju wrote:
> Thanks for the pointer, Vivek. I did check the usual locations, 
> including fe_values.h, but couldn’t locate the exact implementation of 
> the computation for the hessian in the real cell coordinates. (Should 
> involve the product of shape function hessians in reference coordinates, 
> and multiplication with a whole lot of inverse Jacobean terms.)
> 
> Still looking for its exact implementation ….. as I am curious about the 
> computational cost of this implementation.

The various finite element implementations provide the Hessian matrices 
on the reference cell, and these then need to be mapped to the real 
cell. For example, for the usual Q(k) elements, the Hessians on the 
reference cell are pre-computed once at the very beginning (when the 
get_data() function is called by FEValues) and the transformation 
happens in fe_poly.templates.h in the function FE_Poly::fill_fe_values() 
in this piece of code:

   if (flags & update_hessians && cell_similarity != 
CellSimilarity::translation)
     {
       for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
         mapping.transform(make_array_view(fe_data.shape_hessians, k),
                           mapping_covariant_gradient,
                           mapping_internal,
                           make_array_view(output_data.shape_hessians, k));

       for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
         for (unsigned int i = 0; i < quadrature.size(); ++i)
           for (unsigned int j = 0; j < spacedim; ++j)
             output_data.shape_hessians[k][i] -=
               mapping_data.jacobian_pushed_forward_grads[i][j] *
               output_data.shape_gradients[k][i][j];
     }

This accounts for the fact that the Hessian in real space has two parts, 
one that corresponds to the gradient of the Jacobian on the reference 
cell, plus one that is transforming the Jacobian with the gradient of 
the second derivative of the transformation implied by the mapping.

The function transform() is implemented by the individual mappings, but 
in the case of the most common mappings, you'll find it in 
mapping_q_generic.cc -- you just have to pick out the correct one of the 
many versions of the function.

Best
  Wolfgang

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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/be333681-9516-fe3e-867d-8376ede97720%40colostate.edu.

Reply via email to