How to add 2nd iterator to MeshWorker::mesh_loop ?

[Sylvain Mathonnière]

Dear deal.II users,

I am trying to mix step-12 and step-31 for my project.
I have two coupled equations (named eq1 and eq2) that I am solving in parallel just like in step-31 but for one of them I am using the discontinuous Galerkin method like in step-12 (eq1) and for the other one I use standard Galerkin (eq2). I also have two dof_handler like in step-31 one for each equation. It looks like this :

 FESystem<dim>        fe_eq1;        //FE element for equation 1
 FESystem<dim>        fe_eq2;         //FE element for equation 2
 DoFHandler<dim>      dof_handler_eq1;
 DoFHandler<dim>      dof_handler_eq2;

template <int dim> DG_FEM<dim>::DG_FEM ()
: fe_eq1 (FE_DGQ<dim>(1),10), fe_eq2 (FE_Q<dim>(1), 1), dof_handler_eq1 (triangulation), dof_handler_eq2 (triangulation) {}

The problem :

My issue is with the MeshWorker::mesh_loop. Like in step-12, I am looping over cells, faces and boundaries with workers for eq1 like this :

using Iterator = typename DoFHandler<dim>::active_cell_iterator;

auto cell_worker = [&](const Iterator &  cell, ScratchData<dim> &scratch_data, CopyData & copy_data)
  {//my code here needs data from eq2 as well }

MeshWorker::mesh_loop(dof_handler_eq1.begin_active(),
                        dof_handler_eq1.end(),
                        ,/*other arguments*/,
                        cell_worker, /*other workers*/);

However, since my two equations are coupled, I need to access data from eq2 (solution at former time step) but I cannot do it with this iterator here since it is linked only to dof_handler_eq1.

Idealy, I would like to loop using two iterators, one for dof_handler_eq1 and one for dof_handler_eq2.

I know how to do that "traditionnally", I would have something like

auto       eq1_cell       = dof_handler_eq1.begin_active();
const auto end1        = dof_handler_eq1.end();
auto       eq2_cell       = dof_handler_eq2.begin_active();
const auto end2        = dof_handler_eq2.end();

for (;  eq1_cell != end1 &&  eq2_cell != end2 ; ++ eq1_cell , ++ eq2_cell )
{//my code here
}

but with the cell_worker, I do not know.

Question:
How can I modify the cells worker loop (on dof1) initialisation to have a second iterator linked with a second dof_handler looping at the same time ? or how can I add a concurrentely running for loop to the cell_worker with a different iterator ?

I checked in step 16 and 47 as well as in the doc for MeshWorkers and I did not find a solution to this problem. I hope the question is clear enough, do not hesitate to inquires more details if I am missing something important.

Best regards,

Sylvain Mathonnière

[Wolfgang Bangerth]

On 7/12/21 4:18 AM, Sylvain Mathonnière wrote:
>
> _*Question:*_

> How can I modify the cells worker loop (on dof1) initialisation to have a
> second iterator linked with a second dof_handler looping at the same time ? or
> how can I add a concurrentely running for loop to the cell_worker with a
> different iterator ?
>
> I checked in step 16 and 47 as well as in the doc for MeshWorkers and I did
> not find a solution to this problem. I hope the question is clear enough, do
> not hesitate to inquires more details if I am missing something important.

Sylvain,
it's not going to be possible to hack this into MeshWorker::mesh_loop() in an
easy way. But what you can do is this: Inside your assembler, whenever you are
given a cell iterator pointing to a cell in dof_handler_1, you need to find
the corresponding iterator for dof_handler_2. You can do this in the following
way (up to typos and misspelled function names):
DoFHandler::active_cell_iterator cell_2 (cell_1->get_triangulation(),
cell_1->level(),
cell_1->index(),
cell_1->get_dof_handler());

Best
W.

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

[Sylvain Mathonnière]

This sounds like what I need to do indeed. However, I cannot find a function with such arguments in the help of DoFHandler. The closest I could find is in DoFAccessor

DoFAccessor (
const Triangulation< dim, spacedim > *  tria,
const int  level,
const int  index,
const DoFHandler< dim, spacedim > *  dof_handler  )

Is this the one you were refering to ? but then I do not know how to apply it to obtain an iterator "cell_2" on the 2nd cell.

Something like this does not get me an iterator "cell_2" though:

DoFAccessor(cell_1->get_triangulation(),
            cell_1->level(),
            cell_1->index(),
            dof_handler_2)

and In the documentation of the class DoFAccessor, I could not figure out another function to do that.

Maybe I am not using the good class though ?

Best,

Sylvain

[Wolfgang Bangerth]

On 7/13/21 4:02 AM, Sylvain Mathonnière wrote:
> This sounds like what I need to do indeed. However, I cannot find a function

> with such arguments in the help of *DoFHandler*. The closest I could find is
> in *DoFAccessor*
>
> *DoFAccessor (

> const Triangulation< dim, spacedim > *  tria,
> const int  level,
> const int  index,

> const DoFHandler< dim, spacedim > *  dof_handler  )*

>
> Is this the one you were refering to ? but then I do not know how to apply it
> to obtain an iterator "cell_2" on the 2nd cell.

The call was this:

DoFHandler::active_cell_iterator cell_2 (cell_1->get_triangulation(),
cell_1->level(),
cell_1->index(),
cell_1->get_dof_handler());

Note that it's the *iterator* type I'm using here. The class that's behind
this is TriaActiveIterator, which indeed has the constructor you need here.

Best
W>

[Sylvain Mathonnière]

I tried as you suggested but it seems there is a little error still. I found the constructor in TriaActiveIterator (5/8 constructor) though.

I wrote :

typename DoFHandler<dim>::active_cell_iterator cell_2 (cell_1->get_triangulation(), cell_1->level(), cell_1->index(), dof_handler_2);

and the error says (translated from french):

no matching function for the call to « dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> >::TriaActiveIterator(const dealii::Triangulation<2, 2>&, int, int, const dealii::DoFHandler<2, 2>&) »

so this apparently match what I wrote.

A bit further down in the error when enumerating the 8 existing constructor I got

"no known conversion to convert argument 1 from « const dealii::Triangulation<2, 2> » to « const dealii::Triangulation<2, 2>* »

So the constructor is expected a pointer (as the doc of TriaActiveIterator says) but is getting a value. However I checked the doc of get_triangulation() and it  return a reference so I am confused now.

Also for the creation of the second iterator, I have to indicate at some point the dof_handler I want to link to it, so I guess the last argument should be dof_handler_2 rather than cell_1-> get_dof_handler(), otherwise there is no difference between iterator cell_1 and cell_2 as far as I can tell ?

Best,

Sylvain

[Wolfgang Bangerth]

On 7/15/21 2:45 AM, Sylvain Mathonnière wrote:
>
> *typenameDoFHandler<dim>::active_cell_iteratorcell_2 (cell_1->get_triangulation(),
> cell_1->level(), cell_1->index(), dof_handler_2);*

The compiler error tells you that you need pointers, not references. So this
should work:

typename DoFHandler<dim>::active_cell_iterator
cell_2 (&cell_1->get_triangulation(), cell_1->level(), cell_1->index(),
&dof_handler_2);

Best
W.

[Jean-Paul Pelteret]

To add to what Wolfgang has already said, there’s this entry in our FAQ on this topic:

https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#can-i-convert-triangulation-cell-iterators-to-dofhandler-cell-iterators

--
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/53c90809-04b1-c2f0-3d35-ec11726fa46f%40colostate.edu.

[Sylvain Mathonnière]

Adding the & did the trick indeed. Thank you. Thank you as well for the FAQ entry on this. I completely overlooked it...

Best,

Sylvain