Error while implementing FEFieldFunction to calculate the gradient of velocity on the bottom line of backward facing step flow 2D

[Rajeshwari Kamble]

Hi,

I am using the gradient_list function to calculate the velocity gradient in x direction on the bottom wall of backward facing step and trying to output it using the table handler.

But I am getting an exception of points not found in the coarse grid. The points vector and the gradient vector size is 288.

Exception on processing:

--------------------------------------------------------
An error occurred in line <4343> of file </home/kambler/soft/dealii/dealii/source/grid/grid_tools.cc> in function
    std::tuple<std::vector<typename dealii::Triangulation<dim, spacedim>::active_cell_iterator, std::allocator<typename dealii::Triangulation<dim, spacedim>::active_cell_iterator> >, std::vector<std::vector<dealii::Point<dim, double>, std::allocator<dealii::Point<dim, double> > >, std::allocator<std::vector<dealii::Point<dim, double>, std::allocator<dealii::Point<dim, double> > > > >, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > > > dealii::GridTools::compute_point_locations(const dealii::GridTools::Cache<dim, spacedim>&, const std::vector<dealii::Point<spacedim> >&, const typename dealii::Triangulation<dim, spacedim>::active_cell_iterator&) [with int dim = 2; int spacedim = 2; typename dealii::Triangulation<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >; typename dealii::Triangulation<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >]
The violated condition was:
    std::get<3>(cqmp).size() == 0
Additional information:
    The point <0 0> could not be found inside any of the subcells of a coarse grid cell.
--------------------------------------------------------
An error occurred in line <4343> of file </home/kambler/soft/dealii/dealii/source/grid/grid_tools.cc> in function
    std::tuple<std::vector<typename dealii::Triangulation<dim, spacedim>::active_cell_iterator, std::allocator<typename dealii::Triangulation<dim, spacedim>::active_cell_iterator> >, std::vector<std::vector<dealii::Point<dim, double>, std::allocator<dealii::Point<dim, double> > >, std::allocator<std::vector<dealii::Point<dim, double>, std::allocator<dealii::Point<dim, double> > > > >, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > > > dealii::GridTools::compute_point_locations(const dealii::GridTools::Cache<dim, spacedim>&, const std::vector<dealii::Point<spacedim> >&, const typename dealii::Triangulation<dim, spacedim>::active_cell_iterator&) [with int dim = 2; int spacedim = 2; typename dealii::Triangulation<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >; typename dealii::Triangulation<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >]
The violated condition was:
    std::get<3>(cqmp).size() == 0
Additional information:
    The point <17 0> could not be found inside any of the subcells of a coarse grid cell.

PS: I am newer to using deal ii and I checked the discussion on the question   Usage of FEFieldFunction.vector_value_list on a parallel::distributed::Triangulation posted earlier but still not able to
find the mistake.

[Daniel Arndt]

Rajeshwari,

We need a little bit more information to be able to help you.

- Which version of the library are you using?

- What does the mesh you are using look like?

- Do you see this problem also when running in serial or only in parallel?

It would also help to reduce your test case to a size as small as possible if you want to share it with us.

Best,

Daniel

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/73791db1-68b0-4f49-832b-5896b7563876%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Rajeshwari Kamble]

Thank you for your response Professor Daniel.

I am currently using the 9.2.0-pre version of the library and I have attached the gmsh files of the mesh as well as the screenshot.The error occurs while running in parallel but in case of serial computing, it doesn't display anything after setting up the active cells and the degrees of freedom.

Here is the function to compute the gradient :

template<int dim>
void GLSNavierStokesSolver<dim>::reattachment_length()
{
  TimerOutput::Scope t(computing_timer, "reattachment_length");
  QGauss<dim-1>  face_quadrature_formula(degreeQuadrature_+1);
  const MappingQ<dim>      mapping (degreeVelocity_,femParameters.qmapping_all);
  const int n_q_points = face_quadrature_formula.size();
  const FEValuesExtractors::Vector velocities (0);
  const FEValuesExtractors::Scalar pressure (dim);
  std::vector<double> pressure_values(n_q_points);
  std::vector<Tensor< 2, dim> > velocity_gradients(n_q_points);
  std::vector<double> grad;
  Tensor<1,dim> normal_vector;
  Tensor<2,dim> fluid_stress;
  Tensor<2,dim> fluid_pressure;
  Tensor<1,dim> force;
  std::vector< std::vector< Point< dim > > > qpoints;
   std::vector< std::vector< unsigned int > > maps ;
   std::vector< typename DoFHandler<dim>::active_cell_iterator >   cells;
//  double v_grad;
//  double x;
  FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature_formula, update_values
                                    | update_quadrature_points | update_gradients | update_JxW_values
                                    | update_normal_vectors );

   typename DoFHandler<dim>::active_cell_iterator
   cell = dof_handler.begin_active(), endc = dof_handler.end();
   const unsigned int   component = 0 ;
   std::vector< Point< dim > > p(288);
   const typename Triangulation<dim>::active_cell_iterator        cell_hint;

   for (unsigned int k=0;k<288;k+=1)
   {
       p[k]=Point<dim>(k/10,0);
   }
   GridTools::Cache<dim> ob1 (triangulation,mapping);
   GridTools::compute_point_locations(ob1,p,cell_hint);
   for (; cell!=endc; ++cell)
    {
     if (cell->is_locally_owned())
      {
       for (unsigned int face=0; face < GeometryInfo<dim>::faces_per_cell; face++)
        {
         if (cell->face(face)->at_boundary())
          {
            fe_face_values.reinit (cell, face);
            Functions::FEFieldFunction< dim, DoFHandler<dim>, TrilinosWrappers::MPI::Vector> ob(dof_handler,present_solution,mapping);
            ob.compute_point_locations(p,cells,qpoints,maps);
            ob.gradient_list (p,gradients,component);

          }
         }
       }
     }
   {
     std::cout << std::endl;
     TableHandler table;
     table.add_value("u_x",gradients[0][0]);
     table.add_value("u_y",gradients[0][1]);

      std::cout << "+------------------------------------------+"  << std::endl;
      std::cout << "|  Gradient  summary                          |" << std::endl;
      std::cout << "+------------------------------------------+"  << std::endl;
      table.write_text(std::cout);
   }

        gradient_tables.add_value("time",simulationControl.getTime());
        gradient_tables.add_value("u_x",gradients[0][0]);
        gradient_tables.add_value("u_y",gradients[0][1]);

}

[Daniel Arndt]

Rajeshwari,

The purpose of FEFieldFunction is that you don't have to loop through the cells yourself, but that you can just ask for the value at any given point (in global coordinates) as long as the cell

the point is contained in is not artifical. Internally, FEFieldFunction::gradient_list() calls FEFieldFunction::compute_point_locations() which also fails if there is a point contained in artifical cell.

Because of that and since the points you are trying to evaluate your solution at are the same for all processes, it is to be expected that your code throws an exception.

If I were you, I would call GridTools::distributed_compute_point_locations_try_all(https://www.dealii.org/current/doxygen/deal.II/namespaceGridTools.html#a5a845642f4d205352931267b58055d62)

first to get the set of points which are contained in locally owned or ghosted cells and then pass this subset to FEFieldFunction::gradient_list().

Of course, this means that some information GridTools::distributed_compute_point_locations_try_all() provides are not used and recomputed inside FEFieldFunction::gradient_list().

In case this turns to be too costly, you could replicate the implementation of the latter function to use GridTools::distributed_compute_point_locations_try_all() directly instead.

Then, you could also make use of the fact that all the points you are interested in are located at the boundary.

I would definitely try if the first approach is sufficient first, though.

Best,

Daniel