Manually setting Cell dof indices

[Gordan Šegon]

Dear deal.ii team,

I am trying to remap dof indices of certain cells before assembly using the following (sketch code):

    ...
    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
    // Iterate over cells and remap dofs on marked cells
    for (auto& cell : dof_handler.active_cell_iterators()) {
        std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

        if ( /* cell marked for remapping */ ){
            fe_values.reinit(cell);
            cell->get_dof_indices(local_dof_indices);

    for (auto &local_dof_index: local_dof_indices)
local_dof_index = custom_remap_function(local_dof_index);
   
            cell->set_dof_indices(local_dof_index);
        }

    }

Function custom_remap_function implements the user defined mapping.

Compiling produces the following error:

error: passing ‘const dealii::DoFCellAccessor<2, 2, false>’ as ‘this’ argument discards qualifiers [-fpermissive]
  128 |             cell->set_dof_indices(local_dof_indices);

Could you perhaps comment on the outlined strategy and propose a solution - if reasonable. Hopefully I am just missing some trivialities or misusing C++.

At the moment, I implemented the re-mapping during assembly which works fine, but for various reasons would like to perform it before. The motivation for doing this is to (virtually) rotate/move a part of the domain.

~ Gordan Segon

P.S. I am using deal.ii for the past 2 years and would like to thank you all for making the library open source and for providing all the learning materials freely available!

[Wolfgang Bangerth]

On 7/23/23 06:16, Gordan Šegon wrote:
>
> I am trying to remap dof indices of certain cells before assembly using the
> following (sketch code):
>
>     ...
>     const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
>     // Iterate over cells and remap dofs on marked cells
>     for (auto& cell : dof_handler.active_cell_iterators()) {
>         std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
>
>         if ( /* cell marked for remapping */ ){
>             fe_values.reinit(cell);
>             cell->get_dof_indices(local_dof_indices);
>
>     for (auto &local_dof_index: local_dof_indices)
> local_dof_index = custom_remap_function(local_dof_index);
>
>             cell->set_dof_indices(local_dof_index);
>         }
>     }

This will not work unless you have a discontinuous FE because you will be
touching DoFs twice.

A better way is to simply set up a permutation vector mapping old to new DoF
indices and calling DoFHandler::renumber_dofs().

Best
W.

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

[Gordan Šegon]

Hello,

Thank you for the reply. I'm not sure I follow though. 

This will not work unless you have a discontinuous FE because you will be
touching DoFs twice.

What is meant by touching the dofs twice?

The general idea was:

1. The mesh is comprised of stationary and rotary parts. For example, a 2D cylinder is the stationary part and the inner circle of the cylinder is the rotary part. 

2. Vertices on the rotating boundary are spaced equidistant along the boundary.

3. For all cells that have vertices on the rotating boundary, and that are part of the either stationary or rotary part of the mesh, remap the local_dof_index to correspond to a dof_index on the rotating boundary an integer offset away from the local_dof_index when numbered along the periphery of the rotating boundary. This is done during assembly into the global matrix. 

Visualization: https://www.youtube.com/watch?v=Vk__d9avLyE . The visualization of the solution on the first layer of elements on the rotating boundary that are found in the stationary part is in this case incorrect. This is due to no local_dof_index remapping implemented during visualization. Never though, the field is correct in all other cells.

The concept described works but has drawbacks and is hard to debug. One drawback is that the resolution of the rotation is prescribed by the distance between vertices on the rotary boundary and the distances have to be equidistant. Next I will next try to simplify the process by having 2 separated triangulations, performing rotation on one, merging and proceeding as usual. This should work correctly as a reference solution before trying more advanced methods. Eventually, I'd like to have as an input a rotation angle of the rotating mesh, rather than an integer offset. 

A better way is to simply set up a permutation vector mapping old to new DoF
indices and calling DoFHandler::renumber_dofs().

In this context I think this wouldn't work, since this would renumber the dofs for both the outside and inside cells sharing the vertex. If I'm not mistaken.

Is there any standard way of implementing wanted behavior without touching the mesh?

Best regards!

- Gordan

[Wolfgang Bangerth]

> /This will not work unless you have a discontinuous FE because you will be
> touching DoFs twice./

> What is meant by touching the dofs twice?

In your original code snippet, you loop over all cells and on each cell
you renumber all DoF indices of DoFs that live on that cell.

But take a mesh with 4 cells and a Q1 element as an example. Now
consider the following picture:

0--1--2
| | |
3--4--5
| | |
6--7--8

Say you're starting with the bottom left cell. Then you're renumbering
DoFs 3,4,6,7 to whatever new numbers you want to assign them. Say it is
0,1,2,3, and you get the following picture:

0--1--2
| | |
0--1--5
| | |
2--3--8

The next cell is the bottom right cell. Note that if you apply the same
algorithm again, you would renumber the two DoFs on the interface
(originally 4,7, now 1,3) *a second time*.

This is not correct. You have to keep track which DoFs you have already
translated -- which is of course exactly what
DoFHandler::renumber_dofs() does.

> The general idea was:
> 1. The mesh is comprised of stationary and rotary parts. For example, a
> 2D cylinder is the stationary part and the inner circle of the cylinder
> is the rotary part.
> 2. Vertices on the rotating boundary are spaced equidistant along the
> boundary.
> 3. For all cells that have vertices on the rotating boundary, and that
> are part of the either stationary or rotary part of the mesh, remap the
> local_dof_index to correspond to a dof_index on the rotating boundary an

> integer /offset /away from the local_dof_index when numbered along the

> periphery of the rotating boundary. This is done during assembly into
> the global matrix.

I don't understand how that is supposed to work. You assume that there
is an abstract degree of freedom and that cells can arbitrarily
reference that degree of freedom. But that is not the case in deal.II.
For continuous elements such as Q1, the degree of freedom is associated
with a specific vertex, and each of the adjacent cells use this DoF. If
two cells share a vertex, then they also share a unique DoF -- you
cannot break this.

There are ways to work around this. For example, you could use the hp
framework (step-46) and use one field on the inside and one field on the
outside, and then compute constraints that ensure what you want to do.
Or you use discontinuous elements right away. But you can't use a single
field with continuous elements for both inside and outside.


> /A better way is to simply set up a permutation vector mapping old to
> new DoF
> indices and calling DoFHandler::renumber_dofs()./

> In this context I think this wouldn't work, since this would renumber
> the dofs for both the outside and inside cells sharing the vertex. If
> I'm not mistaken.

Yes. But as I mention above, you can't do it separately.

[Gordan Šegon]

In your original code snippet, you loop over all cells and on each cell
you renumber all DoF indices of DoFs that live on that cell.

The loop is over all cells, but only the marked cells would require a modification of the local_dof_indices. My assumption was that each cell holds a copy of local_dof_indices that can be modified without affecting other cells and/or renumbering dofs. By renumbering I mean to change vertex -> dof correspondence.

But take a mesh with 4 cells and a Q1 element as an example. Now
consider the following picture:

0--1--2
| | |
3--4--5
| | |
6--7--8

Say you're starting with the bottom left cell. Then you're renumbering
DoFs 3,4,6,7 to whatever new numbers you want to assign them. Say it is
0,1,2,3, and you get the following picture:

0--1--2
| | |
0--1--5
| | |
2--3--8

The next cell is the bottom right cell. Note that if you apply the same
algorithm again, you would renumber the two DoFs on the interface
(originally 4,7, now 1,3) *a second time*.

This is not correct. You have to keep track which DoFs you have already
translated -- which is of course exactly what
DoFHandler::renumber_dofs() does.

Right, as stated above, this is not wanted. Rather the remapping of certain values in certain marked cell's local_dof_indices is needed. Basically this would distribute the local cell system matrix to different parts of the global system matrix, coupling a given cell to some other cells which are not necessary it's immediate neighbors, depending on the content of local_dof_indices. Under certain assumptions, this should be the same as rotating a portion of the mesh.

I don't understand how that is supposed to work. You assume that there
is an abstract degree of freedom and that cells can arbitrarily
reference that degree of freedom. But that is not the case in deal.II.
For continuous elements such as Q1, the degree of freedom is associated
with a specific vertex, and each of the adjacent cells use this DoF. If
two cells share a vertex, then they also share a unique DoF -- you
cannot break this.

The association of vertex to dof is not wanted to be broken. As mentioned no renumbering is wanted, but only the local_dof_indices modified on per cell basis. However, I see now that perhaps each cell doesn't have it's own copy of local_dof_indices and that changing values local_dof_indices of one cell would affect the local_dof_indices of another?

At the moment the following is done to achieve the "rotation":

For each cell:

2. Assemble local cell_matrix as usual.

3. If cell is marked for "rotation" modify local_dof_indices otherwise do nothing.

4. Distribute the local cell_matrix to the global system_matrix.

And the initial question was if step 3. can be done before local assembly.

I attached 4 pictures that might illustrate the idea easier. The illustration shows how local_dof_indices would be modified for rotation offsets 0, 1, 2. Similarly, local_dof_indices for C2 and C3 would also be modified but not shown in the illustration. C4, C5, C6 are "stationary" and no modification to local_dof_indices is made. 

There are ways to work around this. For example, you could use the hp
framework (step-46) and use one field on the inside and one field on the
outside, and then compute constraints that ensure what you want to do.
Or you use discontinuous elements right away. But you can't use a single
field with continuous elements for both inside and outside.

Thank you for the suggestion, I will definitely look into that. At the moment I am not really happy with the described concept, it seems unnecessarily obscure and "hacky" and probably there is a better way.

Best regards!

- Gordan

[Wolfgang Bangerth]

On 8/21/23 09:36, Gordan Šegon wrote:
> /
> /

> The association of vertex to dof is not wanted to be broken. As mentioned no
> renumbering is wanted, but only the local_dof_indices modified on per cell
> basis. However, I see now that perhaps each cell doesn't have it's own copy of
> local_dof_indices and that changing values local_dof_indices of one cell would
> affect the local_dof_indices of another?

Yes. The DoF (of a continuous finite element) is associated with the vertex,
and the vertex is part of all adjacent cells. You cannot break these associations.