MappingQEulerian produces strange results in 2-d

[jppelteret]

From testing #2205, I've found that using MappingQEularian in conjunction with DataOut in 2-d produces an incompatible displacement field. Attached is an example of this result that illustrates tearing that is not really there. This is a single material with the geometry formed by the GridGenerator::hyper_rectangle function, so all elements are definitely properly connected. Applying the warp function in Paraview renders the right result.

[bangerth]

Interesting. Is the displacement field based on a discontinuous finite element?

[jppelteret]

No it is not - we use FE_Q's. Also this works perfectly fine in 3-d, which is why I hadn't noticed it before :-/

[bangerth]

What's the way to reproduce this problem?

[jppelteret]

I'll see if I can bash out a small testcase later this evening :-)

[jppelteret]

Unfortunately I haven't been able to reproduce the issue thus far. The test below works for polynomial degree 1 and 2. I can't imagine that the pathology is somehow specifically related to the solution generated in step-44. Maybe what I'll try next is to have a multi-component DoFHandler as is the case in that tutorial.

#include <deal.II/base/tensor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>

#include <iostream>
#include <fstream>
#include <vector>
#include <string>

using namespace dealii;

template<int dim>
void MapQE (void)
{
  const unsigned int degree = 2; // Works for degree 1 as well
  const unsigned int n_global_refinements = 2; // Use 3 refinement levels for degree 1

  Triangulation<dim> triangulation;
  GridGenerator::hyper_cube (triangulation,-1,1);
  triangulation.refine_global (n_global_refinements); // Use 3 refinement levels for degree 1

  const FESystem<dim> fe (FE_Q<dim> (degree), dim);
  DoFHandler<dim> dof_handler (triangulation);
  dof_handler.distribute_dofs (fe);


  Vector<double> soln (dof_handler.n_dofs());

  // Apply grid deformation
  Tensor<2,dim> grad_u;
  grad_u[0][0] = 0.1;
  grad_u[0][1] = 0.2;
  grad_u[1][0] = 0.3;
  grad_u[1][1] = 0.4;
  if (dim == 3)
  {
    grad_u[0][2] += 0.5;
    grad_u[2][0] += 0.6;
    grad_u[1][2] += 0.7;
    grad_u[2][1] += 0.8;
    grad_u[2][1] += 0.9;
  }

  Quadrature<dim> q_cell (fe.get_unit_support_points());
  FEValues<dim> fe_values (fe, q_cell, update_q_points);
  typename DoFHandler<dim>::active_cell_iterator
  cell = dof_handler.begin_active(),
  endc = dof_handler.end();
  for (; cell!=endc; ++cell)
  {
    fe_values.reinit(cell);
    std::vector<types::global_dof_index> local_dof_indices(fe.dofs_per_cell);
    cell->get_dof_indices(local_dof_indices);

    for (unsigned int I=0; I<fe.dofs_per_cell; ++I)
    {
      const unsigned int component_I = fe.system_to_component_index(I).first;
      const Point<dim> X = fe_values.get_quadrature_points()[I];
      const Tensor<1,dim> x = grad_u*X;
      soln(local_dof_indices[I]) = x[component_I];
    }
  }

  // Without mapping
  {
    DataOut<dim> data_out;
    data_out.attach_dof_handler (dof_handler);
    data_out.add_data_vector (soln, 
                              std::vector<std::string> (dim, "displacement"),
                              DataOut<dim>::type_dof_data,
                              std::vector<DataComponentInterpretation::DataComponentInterpretation> (dim, DataComponentInterpretation::component_is_part_of_vector));
    data_out.build_patches (degree);
    const std::string filename
     = std::string("solution-")
     + Utilities::int_to_string(dim)
     + std::string("d-001.vtk");
    std::ofstream output (filename);
    data_out.write_vtk (output);
  }

  // With mapping
  {
    DataOut<dim> data_out;
    data_out.attach_dof_handler (dof_handler);
    data_out.add_data_vector (soln, 
                              std::vector<std::string> (dim, "displacement"),
                              DataOut<dim>::type_dof_data,
                              std::vector<DataComponentInterpretation::DataComponentInterpretation> (dim, DataComponentInterpretation::component_is_part_of_vector));
    MappingQEulerian<dim> q_mapping(degree, dof_handler, soln);
    data_out.build_patches(q_mapping, degree);
    const std::string filename
     = std::string("solution-")
     + Utilities::int_to_string(dim)
     + std::string("d-002.vtk");
    std::ofstream output (filename);
    data_out.write_vtk (output);
  }
}

int main()
{
  deallog.depth_console(0);
  std::cout << std::setprecision(10);

  MapQE<2>();
  // MapQE<3>();
}

[bangerth]

Thanks for trying to come up with something that helps us reproduce this in a smaller setting. From the pictures, is it conceivable that in the output the y-coordinate of every point is computed correctly but the x-coordinate is not? Just trying to poke in the dark here...

[jppelteret]

Spurred on by a new question in the forum, here's another attempt at a MWE. It produces a displaced grid for a single and multi-field "problem".

Edit: Removed erroneous code and results. See updated code in this comment.

[davydden]

@jppelteret

// Why does project() not work as expected?

what exactly does not work here? The comment part of the code looks ok to me.

[davydden]

@jppelteret as for the project in three-field case, what you do looks fine as well:

BlockVector tmp(dofs_per_block);
tmp.collect_sizes();
// Set displacement field u
const Displacement u_function (first_u_component,n_components);
VectorTools::project (dof_handler,
constraints,
QGauss(degree+2),
u_function,
tmp);

however it should also be possible just to wrap everything (displacement, pressure, dilatation) into a single Function<dim> and call project once.

[jppelteret]

@davydden Re your last comment: Sure, I agree. I picked this up from some other tests I've been writing, so its a minimal transfer from another test-case to this.

[jppelteret]

@davydden Thanks for doing a sanity check on the code I posted.

Edit: Removed erroneous results. See updated code in this comment.

[davydden]

@jppelteret let's chat about this on Monday unless someone will give a solution during the weekend.

[jppelteret]

Sure, as I said it was a quick attempt to see if I could concretely nail down the issue. If someone can spot an obvious flaw in what I posted I'd be happy to hear about it! But, again, I myself haven't put too much thought into what's going on here.

[jppelteret]

Found the problem:
This

values[first_u_component+0] = X[0]*X[0]*X[1];
values[first_u_component+1] = X[1]*X[1]*X[0];

should be

values[p][first_u_component+0] = X[0]*X[0]*X[1];
values[p][first_u_component+1] = X[1]*X[1]*X[0];

I'll update the MWE accordingly and double check the results.

[jppelteret]

Ok, here's a new piece of code that successfully reproduces the problem. Compare the displaced result solution-one_field-without_mapping-2d-deg-2.vtk to solution-one_field-with_mapping-2d-deg-2.vtk. In the latter case there's tearing near the centre of the displaced geometry.

No Eulerian map (postprocessed with Paraview warp filter), degree=2:

With Eulerian map, degree=2:

Updated MWE / potential unit test:

// Checks that MappingQEularian produces the correct displacement field
// for both single and multi-field systems

// #include "../tests.h"

#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <iostream>
#include <fstream>
#include <cmath>

using namespace dealii;

template <int dim>
class Displacement : public Function<dim>
{
public:
  Displacement (const int first_u_component,
                const int n_components)
    : Function<dim>(n_components),
      first_u_component (first_u_component)
  { }
  virtual void
  vector_value_list (const std::vector< Point<dim> > &points,
                     std::vector< Vector<double> >   &values) const
  {
    for (unsigned int p=0; p<points.size(); ++p)
      {
        const Point<dim> &X = points[p];
        values[p][first_u_component+0] = X[0]*X[0]*X[1];
        values[p][first_u_component+1] = X[1]*X[1]*X[0];
      }
  }

private:
  const int first_u_component;
};

template<int dim>
void
test_one_field_NB (const int degree,
                   const int n_initial_refinements = 2)
{
  Triangulation<dim> triangulation;
  const FESystem<dim> fe(FE_Q<dim>(degree), dim); // displacement
  const QGauss<dim> qf_cell (degree+1);
  DoFHandler<dim> dof_handler(triangulation);

  GridGenerator::hyper_cube(triangulation,0,1);
  triangulation.refine_global(n_initial_refinements);
  dof_handler.distribute_dofs(fe);

  ConstraintMatrix constraints;
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  constraints.close();

  enum
  {
    u_block = 0
  };

  const unsigned int n_blocks = 1;
  const unsigned int n_components = dim;
  const unsigned int first_u_component = 0;

  const FEValuesExtractors::Vector u_fe(first_u_component);

  Vector<double> solution_total(dof_handler.n_dofs());

  {
    Vector<double> tmp(dof_handler.n_dofs());
    // Set displacement field u
    const Displacement<dim> u_function (first_u_component,n_components);
    VectorTools::project (dof_handler,
                          constraints,
                          QGauss<dim>(degree+2),
                          u_function,
                          solution_total);
  }

  // Output solution field
  for (unsigned int output_with_mapping = 0;
       output_with_mapping<=1; ++output_with_mapping)
  {
    DataOut<dim> data_out;
    std::vector<DataComponentInterpretation::DataComponentInterpretation>
    data_component_interpretation(dim,
                                  DataComponentInterpretation::component_is_part_of_vector);
    std::vector<std::string> solution_name(dim, "displacement");
    data_out.attach_dof_handler(dof_handler);
    data_out.add_data_vector(solution_total,
                             solution_name,
                             DataOut<dim>::type_dof_data,
                             data_component_interpretation);

    MappingQEulerian<dim> q_mapping(degree, dof_handler, solution_total);
    if (output_with_mapping == true)
      data_out.build_patches(q_mapping, degree);
    else
      data_out.build_patches(degree);
    std::ostringstream filename;
    filename
      << "solution-one_field-NB-"
      << (output_with_mapping == true ? "with_mapping-" : "without_mapping-")
      << dim << "d-deg-" << degree
      << ".vtk";
    std::ofstream output(filename.str().c_str());
    data_out.write_vtk(output);
  }
}

template<int dim>
void
test_one_field (const int degree,
                const int n_initial_refinements = 2)
{
  Triangulation<dim> triangulation;
  const FESystem<dim> fe(FE_Q<dim>(degree), dim); // displacement
  const QGauss<dim> qf_cell (degree+1);
  DoFHandler<dim> dof_handler(triangulation);

  GridGenerator::hyper_cube(triangulation,0,1);
  triangulation.refine_global(n_initial_refinements);
  dof_handler.distribute_dofs(fe);

  ConstraintMatrix constraints;
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  constraints.close();

  enum
  {
    u_block = 0
  };

  const unsigned int n_blocks = 1;
  const unsigned int n_components = dim;
  const unsigned int first_u_component = 0;

  const FEValuesExtractors::Vector u_fe(first_u_component);

  std::vector<unsigned int> block_component(n_components, u_block); // Displacement
  std::vector<types::global_dof_index> dofs_per_block (n_blocks);
  DoFTools::count_dofs_per_block(dof_handler,
                                 dofs_per_block,
                                 block_component);
  BlockVector<double> solution_total(dofs_per_block);
  solution_total.collect_sizes();

  {
    BlockVector<double> tmp(dofs_per_block);
    tmp.collect_sizes();
    // Set displacement field u
    const Displacement<dim> u_function (first_u_component,n_components);
    VectorTools::project (dof_handler,
                          constraints,
                          QGauss<dim>(degree+2),
                          u_function,
                          tmp);
    solution_total.block(u_block) = tmp.block(u_block);
  }

  // Output solution field
  for (unsigned int output_with_mapping = 0;
       output_with_mapping<=1; ++output_with_mapping)
  {
    DataOut<dim> data_out;
    std::vector<DataComponentInterpretation::DataComponentInterpretation>
    data_component_interpretation(dim,
                                  DataComponentInterpretation::component_is_part_of_vector);
    std::vector<std::string> solution_name(dim, "displacement");
    data_out.attach_dof_handler(dof_handler);
    data_out.add_data_vector(solution_total,
                             solution_name,
                             DataOut<dim>::type_dof_data,
                             data_component_interpretation);

    Vector<double> soln(solution_total.size());
    for (unsigned int i = 0; i < soln.size(); ++i)
      soln(i) = solution_total(i);
    MappingQEulerian<dim> q_mapping(degree, dof_handler, soln);
    if (output_with_mapping == true)
      data_out.build_patches(q_mapping, degree);
    else
      data_out.build_patches(degree);
    std::ostringstream filename;
    filename
      << "solution-one_field-"
      << (output_with_mapping == true ? "with_mapping-" : "without_mapping-")
      << dim << "d-deg-" << degree
      << ".vtk";
    std::ofstream output(filename.str().c_str());
    data_out.write_vtk(output);
  }
}

template<int dim>
void
test_three_field (const int degree,
                  const int n_initial_refinements = 2)
{
  Triangulation<dim> triangulation;
  const FESystem<dim> fe(
    FE_Q<dim>(degree), dim,              // displacement
    FE_DGPMonomial<dim>(degree - 1), 1,  // pressure
    FE_DGPMonomial<dim>(degree - 1), 1); // dilatation
  const QGauss<dim> qf_cell (degree+1);
  DoFHandler<dim> dof_handler(triangulation);

  GridGenerator::hyper_cube(triangulation,0,1);
  triangulation.refine_global(n_initial_refinements);
  dof_handler.distribute_dofs(fe);

  ConstraintMatrix constraints;
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  constraints.close();

  enum
  {
    u_block = 0,
    p_block = 1,
    J_block = 2
  };

  const unsigned int n_blocks = 3;
  const unsigned int n_components = dim + 2;
  const unsigned int first_u_component = 0;
  const unsigned int p_component = dim;
  const unsigned int J_component = dim + 1;

  const FEValuesExtractors::Vector u_fe(first_u_component);
  const FEValuesExtractors::Scalar p_fe(p_component);
  const FEValuesExtractors::Scalar J_fe(J_component);

  std::vector<unsigned int> block_component(n_components, u_block); // Displacement
  block_component[p_component] = p_block; // Pressure
  block_component[J_component] = J_block; // Dilatation
  std::vector<types::global_dof_index> dofs_per_block (n_blocks);
  DoFTools::count_dofs_per_block(dof_handler,
                                 dofs_per_block,
                                 block_component);
  BlockVector<double> solution_total(dofs_per_block);
  solution_total.collect_sizes();

  {
    BlockVector<double> tmp(dofs_per_block);
    tmp.collect_sizes();
    // Set displacement field u
    const Displacement<dim> u_function (first_u_component,n_components);
    VectorTools::project (dof_handler,
                          constraints,
                          QGauss<dim>(degree+2),
                          u_function,
                          tmp);
    solution_total.block(u_block) = tmp.block(u_block);
    // Set constant pressure field p_tilde=5
    const ConstantFunction<dim> p_function (5.0, n_components);
    VectorTools::project (dof_handler,
                          constraints,
                          QGauss<dim>(degree+2),
                          p_function,
                          tmp);
    solution_total.block(p_block) = tmp.block(p_block);
    // Set constant dilatation field J_tilde=1.1
    const ConstantFunction<dim> J_function (1.1, n_components);
    VectorTools::project (dof_handler,
                          constraints,
                          QGauss<dim>(degree+2),
                          J_function,
                          tmp);
    solution_total.block(J_block) = tmp.block(J_block);
  }

  // Output solution field
  for (unsigned int output_with_mapping = 0;
       output_with_mapping<=1; ++output_with_mapping)
  {
    DataOut<dim> data_out;
    std::vector<DataComponentInterpretation::DataComponentInterpretation>
    data_component_interpretation(dim,
                                  DataComponentInterpretation::component_is_part_of_vector);
    data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
    data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);
    std::vector<std::string> solution_name(dim, "displacement");
    solution_name.push_back("pressure");
    solution_name.push_back("dilatation");
    data_out.attach_dof_handler(dof_handler);
    data_out.add_data_vector(solution_total,
                             solution_name,
                             DataOut<dim>::type_dof_data,
                             data_component_interpretation);

    Vector<double> soln(solution_total.size());
    for (unsigned int i = 0; i < soln.size(); ++i)
      soln(i) = solution_total(i);
    MappingQEulerian<dim> q_mapping(degree, dof_handler, soln);
    if (output_with_mapping == true)
      data_out.build_patches(q_mapping, degree);
    else
      data_out.build_patches(degree);
    std::ostringstream filename;
    filename
      << "solution-three_field-"
      << (output_with_mapping == true ? "with_mapping-" : "without_mapping-")
      << dim << "d-deg-" << degree
      << ".vtk";
    std::ofstream output(filename.str().c_str());
    data_out.write_vtk(output);
  }
}

int main (int argc,char **argv)
{
  std::ofstream logfile("output");
  deallog << std::setprecision(3);
  deallog.attach(logfile);
  deallog.threshold_double(1.e-10);

  // Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());

  for (unsigned int deg=1; deg<=2; ++deg)
  {
    deallog.push("deg=" + Utilities::to_string(deg));

    deallog.push("dim=2");
    test_one_field_NB<2>(deg);
    test_one_field<2>(deg);
    test_three_field<2>(deg);
    deallog.pop();

    // deallog.push("dim=3");
    // test_one_field_NB<3>(deg);
    // test_one_field<3>(deg);
    // test_three_field<3>(deg);
    // deallog.pop();

    deallog.pop();
  }

  deallog << "OK" << std::endl;

  return 0;
}

[davydden]

@jppelteret so in the updated MWE we only need test_one_field(), right? I suppose we can say for sure that the problem is in MappingQEulerian which gives those tearing as compared to the post-processed results by Paraview? Everything else seems to work as expected. Does the problem appear for degree=1?

[bangerth]

Could this have something to do with this issue: https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#in-my-graphical-output-the-solution-appears-discontinuous-at-hanging-nodes ?

[davydden]

@bangerth AFAICT there are no hanging nodes here.

[jppelteret]

@davydden For the purpose of this particular issue, yes, we only need test_one_field().
@bangerth Denis is correct: There are no hanging nodes (I recognise that the call to DoFTools::make_hanging_nodes() is unnecessary. It was just "defensive" coding :-) And thanks for correcting my spelling in the post title.

Serendipitously the other test test_three_field() has been useful to highlight another bug! It looks like the DataOut is not outputting the correct results in either scenario, and the root cause is something else entirely!

The issue appears to be there's some expectation that one has called DoFRenumbering::component_wise(dof_handler, block_component); for the multi-field problem.
When one adds this line, the solution is then correct, but the tearing is still present for degree>1 and using MappingQEulerian. Here's the output with some DoF renumbering, degree=3 and using MappingQEulerian:

I've attached a complete set of code and results for anyone who is interested to peruse.
MWE_mapping_q_eulerian.zip

[jppelteret]

@davydden I forgot to add that everything does appear to be working as expected for degree=1.

[jppelteret]

@davydden and I believe that this has something to do with build_patches and its interaction with the mapping. The outputted vertices that are coincident with the actual mesh vertices are in the correct place. Its the additional vertices, placed due to the extra subdivisions specified by the argument to build_patches, which are incorrectly placed (--> discontinuous).

[bangerth]

@jppelteret -- do you think you can whittle this down further to just a few elements (say, 2x2)? That would greatly help debug the issue.

[jppelteret]

Sure, I’ll try my best to do so!

…

On 23 Jan 2017, at 17:18, Wolfgang Bangerth ***@***.***> wrote: @jppelteret <https://github.com/jppelteret> -- do you think you can whittle this down further to just a few elements (say, 2x2)? That would greatly help debug the issue. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2206 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AJ-dwpREAIvR-P98tulAKWrq0VVwjapEks5rVNLAgaJpZM4HcABB>.

[jppelteret]

@bangerth Attached is the smallest MWE that I could produce. The most simplified conditions that I could make reproduce the error are

• Degree = 2
• Global refinement level = 2
• Deformation map: Nonlinear

Strangely, the problem does not manifest itself when the global refinement level is 1 (even if a quadratic FE is used). Also, the deformation map has to be nonlinear (or, at least not the linear map I tried).

[kronbichler]

Checking the code, this will be solved by replacing

      data_out.build_patches(q_mapping, degree);

by

      data_out.build_patches(q_mapping, degree, DataOut<dim>::curved_inner_cells);

After tackling #12423 (this one is a manifestation that we should use the new variant without MappingQ), we should also make curved_inner_cells the default in the mapping in case we build more than one patch or in case the mapping defines preserves_vertex_locations.