Behavior of distribute_dofs of class DoFHandler

[Francesco Costanzo]

Hi,

I am writing to check whether or not the behavior I see in the DoFHandler class is intentional or not. If yes, I am hoping for an explanation so that I can better understand the rationale for the status quo.  If not, I'd like to know if I can hope to see fixes in future versions of the library.  I apologize in advance for the lengthy message.

PREMISE
----------
I am working on a fluid-structure interaction problem using an immersed method to study the motion of a solid immersed in a control volume filled with a fluid.

I am using two independent triangulations, one for the control volume and one for the solid domain.  Correspondingly, I define two dof_handlers.

I want to have a global numbering of all the degrees of freedom.  This allows me, among other things, to assemble a single constraint matrix for the entire problem.  Letting dh_f and dh_s denote the dof_handlers for the control volume and the solid domain, respectively, I first distribute dofs over the control volume with

dh_f.distribute_dofs (fe_f);

and then over the solid domain with

dh_s.distribute_dofs (fe_s,offset);

where offset = dh_f.n_dofs ().

In short, I want the numbering of the dofs over the solid domain to start where the numbering of the dofs over the control volume left off.

QUESTION 1
--------------
Let the number of *new* dofs created and then distributed over the solid be "new_solid_dofs". 

Instinctively (or, perhaps, naively) I would have guessed that new_solid_dofs would have been equal to dh_s.n_dofs().  Instead, it turns out that dh_s.n_dofs() is equal to "offset + new_solid_dofs".  At the same time, if one checks what dofs are "active" on a cell by cell basis, one discovers that only the newly created dofs have actually been distributed over the cells of the triangulation, which is the behavior I was expecting.

** Here is my first question: It the above behavior intentional? If yes, why? i.e., why not have new_solid_dofs = dh_s.n_dofs()?

QUESTION 2
--------------
The fact that dh_s.n_dofs() is not equal to the number of new dofs created has some 
undesirable side effects.  Here is one example ... let's

1) distribute dofs as indicated above, i.e., with dh_s.distribute_dofs (fe_s,offset);
2) then renumber dofs with DoFRenumbering::boost::Cuthill_McKee (dh_s).

This causes the numbering of the dofs on the solid domain to start from 0.  While this is consistent with the the Cuthill-McKee renumbering algorithm, I would have thought, again naively, that a renumbering algorithm would work with the dofs that are actually distributed over the cells, as opposed to the range of numbers going from zero to the number of total dofs registered with the dof_handler.  Be that as it may, the overall effect is that executing Cuthill_McKee (dh_s) nullifies the initial objective of enumerating the new dofs starting with "offset".  A work-around is as follows:

DoFRenumbering::boost::Cuthill_McKee (dh_s, true)

This forces the renumbering of the dofs to be reversed relative to the default behavior, and therefore start with the highest index number, which will be among those newly created.  

** So, here is my second (multi-) question ... is this intentional?  The fact that the renumbering of dofs is not executed using the dofs indices that have actually been distributed over the cells seems insidious to me ... would the reversing of the renumbering always work as a possible work-around, even using other renumbering algorithms?  If one where to change the behavior of the DoFHandler class so that only the newly created dofs are registered with a particular dof_handler, would the behavior of the DoFRenumbering::boost::Cuthill_McKee (dh_s) meet my expectation?

-- 
Best Regards,
Francesco Costanzo
Professor of Engineering Science and Mechanics, and
Professor of Mathematics

Address:
  Center for Neural Engineering
  The Pennsylvania State University
  W-315 Millennium Science Complex
  University Park, PA 16802
  USA
Phone:   (814) 863-2030
Fax:     (814) 865-9974
mailto:cost...@engr.psu.edu
http://www.esm.psu.edu/Faculty/Costanzo/

[Timo Heister]

Dear Francesco,

the parameter offset in distribute_dofs() is a leftover feature that
never got tested or used and I am not surprised that you discovered
these bugs. I am sorry about that. We are probably going to remove or
at least deprecate this option in the next release.
Instead of using offset you have two options:
1. Create a FESystem with your two finite elements and use one DoFHandler.
2. Keep your code as is and distribute with offset=0 with two
DoFHandlers. When you create a block system (vectors, matrices) you
have to hand .block(2), etc. to the corresponding functions because
each DoFHandler starts numbering from 0.

Option 1 is probably easier to implement if you want to assemble a big
block system, see step-20 for example. If you decouple your problems
manually, you might want to go with option 2. I think step-31 and
step-35 show the approach 2.

Thanks,
Timo

> --
> 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
>
>



--
Timo Heister
http://www.math.tamu.edu/~heister/

[Wolfgang Bangerth]

> Option 1 is probably easier to implement if you want to assemble a big
> block system, see step-20 for example. If you decouple your problems
> manually, you might want to go with option 2. I think step-31 and
> step-35 show the approach 2.

And step-21.
W.

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

[Francesco Costanzo]

Dear Timo,

Thank you for your reply and for giving me a heads-up concerning things that will not be supported in the future.

Of the two scenarios you describe I think that the second one is the only one that can work in my case since I have no option but to work with two independent triangulations. Not only are these triangulations non-matching but also, in some sense, moving relative to one another.

Thanks again, Francesco.