DG Matrix-free: getting face iterator form

[Michał Wichrowski]

Dear all,

I'm trying to use matrix-free with discontinuous Galerkin elements for variable coefficient problem.  The face integral require access to value of coefficient. For the cell I've done the following thing:

//Vectors of values, stored in object:

  AlignedVector< VectorizedArray<Number> > mu;

  AlignedVector< VectorizedArray<Number> > rho;

//Filling the vectors:

  const unsigned int n_cells = this->data->n_macro_cells();

  rho.resize  (n_cells);

  mu.resize  (n_cells);

  for (unsigned int cell=0; cell<n_cells; ++cell)

    {

  for(unsigned int element=0; element<this->data->n_components_filled(cell) ; ++element){

  const unsigned int  ID=this->data->get_cell_iterator(cell,element,0 )->material_id();

  rho[cell][element] = rho_map.at(ID) ;

  mu[cell][element]  = miu_map.at(ID) ;

  }

    }

For the faces I am trying to do something like that:

//Vectors of values, stored in object:

  AlignedVector< VectorizedArray<Number> > mu_inner_faces;

//Filling the vectors:

  const unsigned int n_faces=this->data->n_inner_face_batches();

  mu_inner_faces.resize(n_faces);

  for (unsigned int face=0; face<n_faces; ++face){

  const internal::MatrixFreeFunctions::FaceToCellTopology

  < VectorizedArray< Number >::n_array_elements > face2cell_info=this->data->get_face_info(face);

  for(unsigned int element=0;

  element<this->data->n_active_entries_per_face_batch(face) ;

  ++element){

  const unsigned int interior_cell = face2cell_info.cells_interior[element];

    const unsigned int exterior_cell = face2cell_info.cells_exterior[element];

  

   ??????????????????????????????????????

  }

  }

If I understood correctly interior_cell and exterior_cell will be numbers of cell.

 I have no idea how to get cell iterators from that information. I even found that MAtrixFree class contains attribute

std::vector< std::pair< unsigned int, unsigned int > >	cell_level_index

but unfortunately it is private and  there is no  getter for that. 

I will also need to get user_index from each face and create from that another AlignedVector of parameters. This operations will be done only once so I do not care much about efficiency.

Best, 

Michal

[Martin Kronbichler]

Dear Michal,

We do not offer convenience access for all data structures on faces yet, so feel free to suggest what you need.

That said, for the access of rho and mu on an element, we provide the following access function:

https://www.dealii.org/current/doxygen/deal.II/classFEEvaluationBase.html#a442f1ee55e6e80a58e19d53c7184c1a7

This FEEvaluationBase::read_cell_data function expects an AlignedVector<VectorizedArray<Number>> as the argument with as many components as there are cell batches (plus ghosted ones for the MPI case). So in your code, you would fill the rho and mu fields by a cell loop (for MPI, please run it for "n_cells = this->data->n_macro_cells() + this->data->n_ghost_cell_batches()". Then you can access that data both on cells and faces, i.e., with FEEvaluation and FEFaceEvaluation.

The implementation resolves the `interior_cells` and `exterior_cells` fields for you.

Retrieving the `user_index` is a bit more involved as we do not provide access to a face iterator.

You could do it approximately as follows (sitting on a face like in your code below):

const unsigned int interior_cell = face2cell_info.cells_interior[element];

const Triangulation<dim>::cell_iterator my_cell =
  this->data->get_cell_iterator(interior_cell / VectorizedArray<Number>::n_array_elements,
                             interior_cell % VectorizedArray<Number>::n_array_elements);
const Triangulation<dim>::face_iterator face_it = my_cell->face(face2cell_info.interior_face_no);

Can you try if that works? If yes, we should add some access functions to MatrixFree. Let us know if you are interested in helping with a pull request.

Best,
Martin

--
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/b3c0fccd-7a76-4d90-a573-74719d281c37%40googlegroups.com.

[Michał Wichrowski]

Dear Martin

We do not offer convenience access for all data structures on faces yet, so feel free to suggest what you need.

1)   functions similar to get_cell_iterator for getting iterator to interior/exterior cell +  face index:

std::pair<DoFHandler< dim >::cell_iterator,  unstigned int>  get_face_iterator(const unsigned int face_batch_number,

          const unsigned int vector_number, const bool iterior=true,

         const unsigned int fe_component=0)

2) getter for cell_level_index so that the data from get_face_info could be used easily (just in case)

 

That said, for the access of rho and mu on an element, we provide the following access function:

https://www.dealii.org/current/doxygen/deal.II/classFEEvaluationBase.html#a442f1ee55e6e80a58e19d53c7184c1a7

This FEEvaluationBase::read_cell_data function expects an AlignedVector<VectorizedArray<Number>> as the argument with as many components as there are cell batches (plus ghosted ones for the MPI case). So in your code, you would fill the rho and mu fields by a cell loop (for MPI, please run it for "n_cells = this->data->n_macro_cells() + this->data->n_ghost_cell_batches()". Then you can access that data both on cells and faces, i.e., with FEEvaluation and FEFaceEvaluation.

The implementation resolves the `interior_cells` and `exterior_cells` fields for you.

Retrieving the `user_index` is a bit more involved as we do not provide access to a face iterator.

You could do it approximately as follows (sitting on a face like in your code below):

const unsigned int interior_cell = face2cell_info.cells_interior[element];
const Triangulation<dim>::cell_iterator my_cell =
  this->data->get_cell_iterator(interior_cell / VectorizedArray<Number>::n_array_elements,
                             interior_cell % VectorizedArray<Number>::n_array_elements);
const Triangulation<dim>::face_iterator face_it = my_cell->face(face2cell_info.interior_face_no);

Can you try if that works? If yes, we should add some access functions to MatrixFree. Let us know if you are interested in helping with a pull request.

 

Ok, it works (at least on active level). I tested that by comparing location of quadrature points (code below). But for implementing get_face_iterator I  think that it would be safer to do it using  cell_level_index (eg. numbering of cell may change )

FEFaceEvaluation<dim, degree_u, degree_u+1+0 , dim, Number> phi_inner(*(this->data),

                                                                       true);

FEFaceValues<dim> fe_face_values(this->data->get_dof_handler().get_fe(),

this->data->get_face_quadrature(),

UpdateFlags::update_quadrature_points  );

  const unsigned int n_faces=this->data->n_inner_face_batches();

  mu_inner_faces.resize(n_faces);

  for (unsigned int face=0; face<n_faces; ++face){

  phi_inner.reinit(face);

  const internal::MatrixFreeFunctions::FaceToCellTopology

  < VectorizedArray< Number >::n_array_elements > face2cell_info=this->data->get_face_info(face);

  for(unsigned int element=0;

  element<this->data->n_active_entries_per_face_batch(face) ;

  ++element){

  const unsigned int interior_cell = face2cell_info.cells_interior[element];

//     const unsigned int exterior_cell = face2cell_info.cells_exterior[element];

  const typename Triangulation<dim>::cell_iterator my_cell =

    this->data->get_cell_iterator(interior_cell / VectorizedArray<Number>::n_array_elements,

                               interior_cell % VectorizedArray<Number>::n_array_elements);

  const typename Triangulation<dim>::face_iterator face_it = my_cell->face(face2cell_info.interior_face_no);

  fe_face_values.reinit(my_cell, face2cell_info.interior_face_no);

  std::vector<Point< dim> > iterator_points = fe_face_values.get_quadrature_points ();

  for(unsigned int point_number=1; point_number<iterator_points.size(); ++point_number){

  Point< dim, VectorizedArray< Number > > point_face =phi_inner.quadrature_point (point_number);

double difference =0;

for(unsigned int d=0; d<dim; ++d){

difference+= std::fabs(iterator_points[point_number](d)-point_face(d)[element]);

}

if (difference>1e-8)

std::cout<<"Diffrence:  " <<difference<< " in quadrature points on faces!"<<std::endl;

  }

  mu_inner_faces[face][element]=1;

  }

  }

 

[Martin Kronbichler]

Dear Michal,

Great. I've opened issue https://github.com/dealii/dealii/issues/8460
for that to keep track of it. We won't get at it immediately but I hope
to do so in the coming weeks.

Best,
Martin