different shape_gradients for MappingQ and MappingQEulerian

[thomas stephens]

I am having trouble understanding how to use MappingQEulerian - I have constructed a codimension-1 sphere by attaching a SphericalManifold to a hyper_sphere triangulation.

So far so good.

Now, I have some surface derivatives I have been computing within the function below called assemble_system()

assemble_system ()
{
  /*{{{*/
 
  /*
   *    [   M      L-2*hL+0.5*d ][ V_n+1 ]   [  0  ]
   *    |                       ||       | = |     |
   *    [ -zn*L          M      ][ H_n+1 ]   [ rhs ]
   *
   *    system_matrix*VH = system_rhs
   *   
   *    system_matrix has size 2*n_dofs x 2*n_dofs,
   *    each block has size n_dofs x n_dofs, and
   *    rhs has size n_dofs
   *
   */
 
 
  MappingQ<dim,spacedim> mapping(2);
 
  system_matrix = 0;
  VH = 0;
  system_rhs = 0;

  const QGauss<dim>  quadrature_formula (2*fe.degree);
  FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
                                    update_values              |   
                                    update_normal_vectors      |   
                                    update_gradients           |   
                                    update_quadrature_points   |
                                    update_JxW_values);

and life has been pretty good.  The subsequent derivatives look like this:

  const FEValuesExtractors::Vector W (0);
  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 
  for (typename DoFHandler<dim,spacedim>::active_cell_iterator
       cell = dof_handler.begin_active(), 
       endc = dof_handler.end();
       cell!=endc; ++cell)
  {
    local_M   = 0;
    local_L   = 0;
    local_hL  = 0;
    local_d   = 0;
    local_rhs = 0;

    fe_values.reinit (cell);

    for (unsigned int i=0; i<dofs_per_cell; ++i)
      for (unsigned int j=0; j<dofs_per_cell; ++j)
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
        {

          local_M(i,j)  += fe_values[W].value(i,q_point)*
                           fe_values[W].value(j,q_point)*
                           fe_values.JxW(q_point);

          local_L(i,j)  += scalar_product(fe_values[W].gradient(i,q_point),
                                          fe_values[W].gradient(j,q_point)
                                         )* fe_values.JxW(q_point);

          local_hL(i,j) += scalar_product(2.0*identity_on_manifold.shape_grad(fe_values.normal_vector(q_point))*
                                          fe_values[W].gradient(i,q_point),
                                          fe_values[W].gradient(j,q_point)
                                         )* fe_values.JxW(q_point);

But now I want to move my mesh, and my strategy has been to compute a deformation vector (of size dof_handler.n_dofs()), and create a new MappingQEulerian at each time step using this updated deformation vector. 

My problem is that I cannot really get off the ground with this strategy because (it seems to me that) once I replace MappingQ with MappingQEulerian, these derivatives are no longer correct.  So, for the modified assemble_system() code I have:

assemble_system ()
{
  /*{{{*/
 
  /*
   *    [   M      L-2*hL+0.5*d ][ V_n+1 ]   [  0  ]
   *    |                       ||       | = |     |
   *    [ -zn*L          M      ][ H_n+1 ]   [ rhs ]
   *
   *    system_matrix*VH = system_rhs
   *   
   *    system_matrix has size 2*n_dofs x 2*n_dofs,
   *    each block has size n_dofs x n_dofs, and
   *    rhs has size n_dofs
   *
   */
 
  /* build a MappingQEulerian that approximates
   * grid_transform : sphere --> ellipsoid. This will be updated by the
   * deformation induced by the velocity field computed at each time step.
  */
 
  /* create an interpolation of the transformation that takes the sphere to the
   * ellipsoid, this will be updated at each mesh movement */
 
  deformation = 0;
  MappingQEulerian<dim,Vector<double>,spacedim> mapping(2, dof_handler, deformation);

 
  system_matrix = 0;
  VH = 0;
  system_rhs = 0;

  const QGauss<dim>  quadrature_formula (2*fe.degree);
  FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
                                    update_values              |   
                                    update_normal_vectors      |   
                                    update_gradients           |   
                                    update_quadrature_points   |
                                    update_JxW_values);

The difference in output looks like this:


From MappingQ:























From MappingQEulerian:






Something is very wrong here - I specifically set deformation=0; before building MappingQEulerian in order to test the null deformation (i.e. the identity), but I still get some silly numbers for H2 (which is computed directly from those derivatives.

What am I missing?  Let me know if there is more information I can provide here.


Thank you,
Tom

[thomas stephens]

Here's one way around this:

At the end of the day, the problem I had is that I did not understand how a Manifold object interacts with a Triangulation and a Mapping object.  I probably still don't, so take this with a grain of salt:

I was creating a Manifold object and then attaching it to a Triangulation, like so  

 static Ellipsoid<dim,spacedim> ellipsoid(a,b,c,center);

 

static SphericalManifold<dim,spacedim> sphere();
GridGenerator::hyper_sphere(triangulation,center, 1);  
triangulation.set_all_manifold_ids(0);
triangulation.set_manifold (0, sphere);

After setting up the DoFHandler and FESystem, a call to 

FEValues<dim,spacedim> fe_values(mapping, fe, quadrature_formula,...);

maps all points from the triangulation -- and quadrature points -- according to the push_forward function in SphericalManifold.  The result is that we can easily compute shape gradients on this surface.

Alternatively, when the Mapping is not MappingQ, but MappingQEulerian, the mapping waits for specific instructions on how to move the quadrature points to the desired shape - that sits in the euler_vector. I was passing a zero euler_vector because I figured that my manifold and triangulation would have figured out what to do with all relevant points.  

The fix here is to get rid of the SphericalManifold altogether, and just start with a triangulation. Then create an euler_vector by interpolating a function which can map arbitrary points to the desired geometry (this is essentially the push_forward function of SphericalManifold if you want a sphere).  This euler_vector will not be zero!  

The whole thing looks like this:

GridGenerator::hyper_sphere(triangulation,center,1.0);

euler_vector.reinit (dof_handler.n_dofs());

VectorTools::interpolate(dof_handler,

                         MapToSphere<spacedim>(r),

                         euler_vector);

MappingQEulerian<dim,Vector<double>,spacedim> mapping(2, dof_handler, euler_vector);

FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula, ...)

Now one can obtain a sphere, compute shape derivatives on it, and then update the euler_vector in order to deform the sphere into something else. 

I suggest to myself to reread the Interplay of UpdateFlags, Mapping, and FiniteElement in FEValues: https://dealii.org/8.4.1/doxygen/deal.II/group__UpdateFlags.html, maybe this will all become clear one day ;)