Extracting a component from Block FESystem and DofHandler

[Spencer Patty]

Dear All,

I have a class which represents a velocity field (I am using a variant of the pressure correction navier stokes model with two phases) and a class which represents a transport model along a provided velocity field.   I am adding in some terms into the velocity model (currently surface tension but soon to be added is the canham-helfrich or willmore energy) accounting for the physics and force balance on the boundary between the two phases.  The approach that we are using is to create a block structure inside the velocity class to incorporate these other physics into the model.  We have a splitting of velocity separate from the pressure terms and have separate dofhandlers for these two groups but now we are lumping all the components that need to be simultaneously solved for with the velocity into the velocity dofhandler and fesystem.  For instance we have the block components in the velocity fesystem:

dim components of velocity,  normal velocity, curvature and willmore energy

so that our velocity fesystem looks like

FESystem<dim> (

                             FESystem<dim>(FE_Q<dim> (parameters.degree_of_fe_velocity),dim), 1,  // velocity

                             FE_Q<dim>(parameters.degree_of_fe_velocity), 1, // normal velocity

                             FE_Q<dim>(parameters.degree_of_fe_velocity), 1, // curvature

                             FE_Q<dim>(parameters.degree_of_fe_velocity), 1 // willmore force

                             )

I have figured out how to construct the block sparsity patterns and matrices etc and solve for this, but at the end of the day this class represents a velocity not a velocity plus other stuff.  This class interfaces with the other class representing transport and I need to be able to provide a velocity, dofhandler and fesystem in such a way that I don't break encapsulation ie the other class doesn't need to know that internally there was a block system and other components... 

The other class expects to be able to query a TrilinosWrappers::MPI::Vector,  FESystem<dim> and DoFHandler<dim>  representing the velocity finite element solution.

Supposing internally to the velocity class I have

solution, a TrilinosWrapper::MPI::BlockVector

fe_system, the FESystem<dim> described above and 

dof_handler, the associated DoFHandler<dim> which is sorted DoFRenumbering::hierarchical and DoFRenumbering::component_wise.

Supposing velocity_block represents the block number associated to the velocity components, I noticed that solution.block{velocity_block} would be a TrilinosWrapper::MPI::Vector of the velocity as desired and it is possible to extract the sub FiniteElement fe_system.base_element{velocity_block} corresponding only to the velocity, but the dof_handler has no such existing capability and if I pass it as is with the extracted FiniteElement into an FEValues, an assertion is triggered when we call fe_values.reinit(cell) saying that the given finite element does not match the dofhandler's finite element...  I don't want to pass the entire system since that requires the other class to know that it only needs the first block components etc.

My question is this:  What would be a best approach to dealing with this?  Is there currently a way to extract a "diagonal block" from the system and reduce the solution, fesystem and dofhandler to only represent that part?  Do I need to create a new FESystem, DoFHandler and vector to provide the functionality I need ?

I hope this is clear enough and enough(and not too much) background to understand what I am asking.  Thanks for any help!

Spencer Patty

[Bruno Turcksin]

Spencer,

I think that passing around the entire system is the way to go but what I would do is to use component mask (http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossComponentMask) and block mask (http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossBlockMask) to work on a subset of the whole system. Instead of having a bunch classes working on small blocks that will be added together at the ends, you have an entire system but each class works on a small block without knowing what the others do. For example, the velocity class can have access to the global dof_handler and a FEValuesExtractors for the velocity. It doesn't need to know what is its block number or how many blocks there are.

Best,
Bruno

[Spencer Patty]

Bruno,

Thanks for responding!  

On Tuesday, August 16, 2016 at 11:16:41 AM UTC-5, Bruno Turcksin wrote:

Spencer,

I think that passing around the entire system is the way to go but what I would do is to use component mask (http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossComponentMask) and block mask (http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossBlockMask) to work on a subset of the whole system. Instead of having a bunch classes working on small blocks that will be added together at the ends, you have an entire system but each class works on a small block without knowing what the others do. For example, the velocity class can have access to the global dof_handler and a FEValuesExtractors for the velocity. It doesn't need to know what is its block number or how many blocks there are.

 

Are there examples of how to do what you describe?  Looking at the documentation, it is not apparent how to do this.

 

Best,
Bruno

So supposing I have the public interface of my velocity class (specified by an abstract velocity interface class from which this velocity model derives) having the following member functions:

      TrilinosWrappers::MPI::Vector & get_velocity()

      FESystem<dim> &  get_fe_velocity()

      DoFHandler<dim> & get_dof_handler_velocity()

how would I go about using the Block or Component Mask to return the proper things?

For velocity, I can return  solution.block{velocity_block}  which gives a TrilinosWrappers::MPI::Vector and I could possibly change the interface to return an FiniteElement<dim> instead of FESystem<dim>  so that I could return fe_system.base_element(velocity_block)  but the mask just allows me to extract dofs into a vector,  how do I get that back into a dofhandler?

My velocity model is created using the factory method which keep to the abstract velocity interface class and I have several other velocity models that all obey this interface above.  The other transport class expect this interface, so it would be very cumbersome to change it for this specific block structure since the others do not use it.  Internally it has a block structure but externally, it represents only a velocity which is one of the blocks.  

[Spencer Patty]

On Tuesday, August 16, 2016 at 12:35:05 PM UTC-5, Spencer Patty wrote:

Bruno,

Thanks for responding!  

On Tuesday, August 16, 2016 at 11:16:41 AM UTC-5, Bruno Turcksin wrote:

Spencer,

I think that passing around the entire system is the way to go but what I would do is to use component mask (http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossComponentMask) and block mask (http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossBlockMask) to work on a subset of the whole system. Instead of having a bunch classes working on small blocks that will be added together at the ends, you have an entire system but each class works on a small block without knowing what the others do. For example, the velocity class can have access to the global dof_handler and a FEValuesExtractors for the velocity. It doesn't need to know what is its block number or how many blocks there are.

 

Are there examples of how to do what you describe?  Looking at the documentation, it is not apparent how to do this.

 

I should say,  I am already using the FEValuesExtractors internally to assemble and solve the system although I never used the Block or Component Mask internally,  my question pertains more to what comes after I have the solution and I want to use it elsewhere.  

[Bruno Turcksin]

Sorry, I misunderstood what you were trying to do. Like you found out,
the problem is that you cannot "extract" the DoFHandler because the
DoFHandler has no block structure and it will have to be rebuilt.
Everything else can be extracted. If your interface returns the
Triangulation and the FESystem or the FiniteElement, you can easily
rebuild the DoFHandler. The only thing that you need to be careful
about is to make sure that you are using the same ordering everywhere.

Best,

Bruno

> --
> 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 a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/u0Q4ByjEjCg/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

On 08/16/2016 09:34 AM, Spencer Patty wrote:
>
> Supposing velocity_block represents the block number associated to the
> velocity components, I noticed that solution.block{velocity_block} would be a
> TrilinosWrapper::MPI::Vector of the velocity as desired and it is possible to
> extract the sub FiniteElement fe_system.base_element{velocity_block}
> corresponding only to the velocity, but the dof_handler has no such existing
> capability and if I pass it as is with the extracted FiniteElement into an
> FEValues, an assertion is triggered when we call fe_values.reinit(cell) saying
> that the given finite element does not match the dofhandler's finite
> element...

Correct -- there is no way to subset a DoFHandler by vector components. That's
because a FESystem is really stored as a collection of simpler elements,
whereas a DoFHandler is not.

> My question is this: What would be a best approach to dealing with this? Is
> there currently a way to extract a "diagonal block" from the system and reduce
> the solution, fesystem and dofhandler to only represent that part? Do I need
> to create a new FESystem, DoFHandler and vector to provide the functionality I
> need ?

The first two yes. If you do it right, i.e., ensure that the velocity
DoFHandler numbers degrees of freedom in the same order as the system
DoFHandler does, then you won't need a separate solution vector -- you can
just use a block of your block vector.

The alternative -- and I think that is what Bruno was suggesting in his first
email -- is to keep the one system DoFHandler and vector. In addition to
those, you would then have to pass something to all consumers of the velocity
field that identifies which components of the FESystem and of the block vector
they need to know about.

Cheers
W.

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

[Corbin Foucart]

Hi all,

I'm doing something very similar to Spencer. The answers provided by Prof. Bangerth and Bruno make sense, up to one question:

The first two yes. If you do it right, i.e., ensure that the velocity
DoFHandler numbers degrees of freedom in the same order as the system
DoFHandler does, then you won't need a separate solution vector -- you can
just use a block of your block vector.

 How would one go about ensuring that the velocity DoFHandler numbers the DOF in the same order as the system DoFHandler does? Is there a way to build the former from the latter?

Best,

Corbin

[Wolfgang Bangerth]

Good question. The answer is that if you just call
DoFHandler::distribute_dofs(), then the order in which the velocities are
numbered in the two DoFHandlers is the same. This order is preserved if you
then call DoFRenumbering::component_wise on one or both of them, and that
guarantees that the elements of vectors match.

Best