Stabilized FEM implementation with bubble function

[Lixing Zhu]

Dear all,

I am trying to implement a stabiliazed weak form (e.g. advection-diffusion) where the stabilization tensor is computed element-wise through a standard bubble: \Pi(1-x_i^2). It seems that FE_Q_Bubbles should provides all I need, but here are two things I am not quite clear about,

1. Is the definition of bubble function in FE_Q_Bubbles in [0,1] span or [-1,1] span? Does it has the standard bubble shape (1 at element center and vanishes at the edges)? The (2x_j-1)^{p-1} part is a little bit confusing to me. The bubble function corresponds to linear element in FE_Q_Bubbles is the standard bubble that I desire, but I am really confusing at the shape of this expression at higher orders.

2. I tried to utilize this class (i.e., FE_Q_Bubbles). One issue is that it takes the virtual nodes of the bubble function into account of the total DOFs, which is not the way we prefer in the stabilization method. I hope the bubble function only provides a local support to estimate the stabilization parameter but not involved into the global linear system. How can I deactive the dofs of bubble function when I distributed the dofs like

dof_handler.distribute_dofs(fe) in the tutorial.

I am pretty new with dealii. Any clue or reference paper is fine.

Thanks you guys in advance.

[Wolfgang Bangerth]

Lixing,

> I am trying to implement a stabiliazed weak form (e.g. advection-diffusion)
> where the stabilization tensor is computed element-wise through a standard
> bubble: \Pi(1-x_i^2). It seems that FE_Q_Bubbles should provides all I need,
> but here are two things I am not quite clear about,
>
> 1. Is the definition of bubble function in FE_Q_Bubbles in [0,1] span or
> [-1,1] span?

Everything in deal.II is on the reference cell define as [0,1]^d.

> Does it has the standard bubble shape (1 at element center and
> vanishes at the edges)? The (2x_j-1)^{p-1} part is a little bit confusing to
> me. The bubble function corresponds to linear element in FE_Q_Bubbles is the
> standard bubble that I desire, but I am really confusing at the shape of this
> expression at higher orders.

For higher orders, you end up with multiple bubble functions, one for each
j=0..dim-1. This wasn't quite clear from the documentation, and I'll submit a
patch later for that.

> 2. I tried to utilize this class (i.e., FE_Q_Bubbles). One issue is that it
> takes the virtual nodes of the bubble function into account of the total DOFs,
> which is not the way we prefer in the stabilization method.

But the behavior is correct -- the class describes a finite element *space*
and that space contains the bubble function. I think that what you are trying
to do is to do static elimination of that degree of freedom right away, and
that can be implemented as well but is not the philosophical view we generally
take in deal.II if you select such an element.

The question is what you want the finite element space to be (i) locally, on
every cell, and (ii) globally. There are a number of tutorials that discuss
these sorts of questions. I would encourage you to read through step-61 and
step-51, for example.

Best
WB

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

[Lixing Zhu]

Many thanks for such an instant response. I'll adopt these ideas and try to eliminate those DOFs related to bubble support following the workflow in the tutorial.

Regards,

Lixing

[Lixing Zhu]

Dear Wolgang,

1. I looked into the step-51 of the tutorial. It does illustrate a paradigm of segregating local DOFs and global DOFs. If I utilize this paradigm, the workflow would solve the local DOFs first (virtual node of the bubble function), which is a block matrix since bubble support from adjacent cells is irrelevant. Then I substitute the solution to the system related to global DOFs (vertices of the standard Lagrange shape function), namely, stabilizing the system. However, this would require certain FEM space in the barycenter of the element (like a FE_Center, in the same fashion as FE_Face), which seems not provided in deal.ii.

2. On the other hand, I am more inclined to use element-wise static condensation (I guess this is somewhat not a rational decision). Is there any example of how to eliminate the virtual DOFs of bubbles from the global system once we have eliminated them locally? I guess this route can reduce the time of implementation since I can at least utilize the FE_Q_Bubbles with its linear realization.

3. The bubbles in FE_Q_Bubbles for a quadratic and higher-order element are not the standard bubbles. Are there any specific references to these bubbles? or how hard it is to modify the bubble function in FE_Q_Bubbles for my own choice?

Thank you for your time in advance.

Best

Lixing

On Tuesday, January 5, 2021 at 1:21:56 PM UTC+8 Wolfgang Bangerth wrote:

[Jean-Paul Pelteret]

Hi Lixing,

Another tutorial that might be of interest to you to look at is step-44. In that tutorial, two discontinuous fields are condensed out in one of two ways (there's a switch to choose which method is applied). The first approach uses the "local condensation" technique, where one block of the global block-system matrix is augmented when doing the condensation. The second approach incorporates condensation into the solver strategy for the linear system by adopting a (nested) global Schur complement; unlike the first technique, this one does not change the values in the stored in the system matrix. In either case, just one field of a three field problem is solved for, and the others computed as a post-processing step.

I think that the first approach might cover the details of one possible way of doing what you list as point (2) "element-wise static condensation". However, the second method (using LinearOperators) is, in my opinion, much easier to implement (although, in the end, both approaches have their advantages and disadvantages).

I hope that this helps a little.

Best,
Jean-Paul

--
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/ba1bb5e4-4344-4f19-b630-8dee8d86dc3dn%40googlegroups.com.

[Wolfgang Bangerth]

> 1. I looked into the step-51 of the tutorial. It does illustrate a paradigm of
> segregating local DOFs and global DOFs. If I utilize this paradigm, the
> workflow would solve the local DOFs first (virtual node of the bubble
> function), which is a block matrix since bubble support from adjacent cells is
> irrelevant. Then I substitute the solution to the system related to global
> DOFs (vertices of the standard Lagrange shape function), namely, stabilizing
> the system. However, this would require certain FEM space in the barycenter of
> the element (like a FE_Center, in the same fashion as FE_Face), which seems
> not provided in deal.ii.

You can do as you suggest, but you can also do the local elimination as part
of the assembly process. For example, assume that there are n regular and m
bubble functions, then you would build a (n+m)x(n+m) local matrix. You can
then eliminate the m bubble functions from this local system because you know
that they do not couple with any other DoF on other cells. You'd do this by
adding multiples of these m lines to the first n lines so that you end up with
a block structure of the local matrix that looks like this

[ A_nn 0]
[ A_mn X]

You can think of this in the following way: Let's assume you start with the matrix

[ B_nn B_nm ]
[ B_mn B_mm ]

Then you need to multiply the second block row by -B_{nm}B_{mm}^{-1} and add
the result to the first block row. You then have
A_nn = B_nn - B_{nm}B_{mm}^{-1}B_{mn}

You'd have to do the same thing to the right hand side vectors, of course.

If you do this, then the global matrix with also have this block triangular
structure, and you can solve the problem for the "regular" DoFs using

B_nn U_n = ...

You *could* also solve for the bubble DoFs using

B_mm U_m = rhs_m - A_mn U_n

but you don't actually have to do that if you don't care about these DoFs.

> 2. On the other hand, I am more inclined to use element-wise static
> condensation (I guess this is somewhat not a rational decision). Is there any
> example of how to eliminate the virtual DOFs of bubbles from the global system
> once we have eliminated them locally?

Not that I know of, but the method above should work.

> 3. The bubbles in FE_Q_Bubbles for a quadratic and higher-order element are
> not the standard bubbles. Are there any specific references to these bubbles?
> or how hard it is to modify the bubble function in FE_Q_Bubbles for my own choice?

What is "standard"? The implementation uses one way to define bubbles, but I'm
sure there are others. None of these would be very difficult to implement once
you understand how the FE_Q_Bubbles class is implemented.

Best
W.

[Lixing Zhu]

Dear W.

Many thanks for your suggestion. I guess I have to implement a new class for my own bubbles at some point.

A follow-up question on (1).

I implemented the local elimination as you suggested and it works. However, this ends up with many extra "redundant" DOFs in the global system. For example, if I am using N*N linear elements for a scalar 2D problem with structured mesh, instead of (N+1)^2 global DOFs, I end up with (N+1)^2 + N^2 global DOFs. The extra DOFs are coming from the bubble supports. Of course, they won't affect the solution at vertices due to the nature of bubble functions, but such an increase of DOFs may require more memory in the solution process. How should I also reduce the size of the global system?

Regards,

Lixing

[Lixing Zhu]

Dear Jean-Paul,

Thanks for your suggestion. I'll try to test with these two approaches.

Regards,

Lixing

[Wolfgang Bangerth]

On 1/7/21 6:58 PM, Lixing Zhu wrote:
>
> I implemented the local elimination as you suggested and it works. However,
> this ends up with many extra "redundant" DOFs in the global system. For
> example, if I am using N*N linear elements for a scalar 2D problem with
> structured mesh, instead of (N+1)^2 global DOFs, I end up with (N+1)^2 + N^2
> global DOFs. The extra DOFs are coming from the bubble supports. Of course,
> they won't affect the solution at vertices due to the nature of bubble
> functions, but such an increase of DOFs may require more memory in the
> solution process. How should I also reduce the size of the global system?

Do you need to? Is memory consumption an issue for you? Unless you have a
concrete problem, it's probably not worth optimizing for.

[Lixing Zhu]

Hi W and Jean-Paul,

Thanks again for your suggestions.

I realized that the bubble functions that I prefer are actually the shape functions from other standard elements. Thus I only need to define two sets of FE (i.g. FE_1 and FE_b). Then I use the iterator from triangulation as

for (const auto &cell : triangulation.active_cell_iterators())

    {

      fe_values_1.reinit(cell); // standard Lagrange elements

      fe_values_b.reinit(cell); // bubble functions

This way I can freely use bubble functions and only construct local K and F of FE_1. However, one issue is that, apparently, triangulation.cell does not support cell->get_dof_indices(local_dof_indices). Is there any way to call the DOFhandler by the cell in triangulation and then access the local_dof_indices?

Regards,

Lixing

[Wolfgang Bangerth]

On 1/7/21 10:22 PM, Lixing Zhu wrote:
>
> You are right. I am not bounded by the memory at this moment. The only
> potential issue I can think of at this stage is the computation of the error
> norm. Technically, if I feed the slots in global RHS corresponding to bubble
> support as zero, the solution at the bubble support should zero. Then the
> error norm would be the same as compared to use solutions at vertices only. Is
> this presumption valid here?

No. You are computing a solution in the space that includes the bubbles. (The
fact that you want to employ a clever way to solve the linear problem does not
change this.) So you have to compute the error with the solution including the
bubble functions, unless you only evaluate the solution at points at which the
bubble functions are zero.

[Wolfgang Bangerth]

On 1/8/21 12:15 PM, Lixing Zhu wrote:
>
> This way I can freely use bubble functions and only construct local K and F of
> FE_1. However, one issue is that, apparently, triangulation.cell does not
> support cell->get_dof_indices(local_dof_indices). Is there any way to call the
> DOFhandler by the cell in triangulation and then access the local_dof_indices?

You want a statement like this (from step-43):

auto cell = darcy_dof_handler.begin_active();
const auto endc = darcy_dof_handler.end();
auto saturation_cell = saturation_dof_handler.begin_active();

for (; cell != endc; ++cell, ++saturation_cell)

Here, we're walking iterators over two DoFHandler objects in parallel. In your
case, one will be the DoFHandler iterator and the other one the Triangulation
iterator. Or, you could do
fe_values_b.reinit
(static_cast<typename Triangulation<dim>::active_cell_iterator>
(cell));
to cast the DoFHandler iterator to a Triangulation iterator.

[Lixing Zhu]

Hi Wolfgang,

Many thanks! The issue is resolved and the code is working.

Regards,

Lixing