How does the step-20 work on spaces of Pn and Pn-1 elements?

[jieli...@gmail.com]

Dear all,

In the step-20, the mixed finite element method is used to solve a 2D Laplace problem, where RT and DG elements are used for the velocity and the pressure respectively. How does it work on spaces of Pn and Pn-1 elements?

Since the velocity is a vector, the main problem becomes representing the two components of the velocity both by Pn elements. Thanks for your help!

Best,

Jie Liu

[Wolfgang Bangerth]

On 09/15/2017 09:51 AM, jieli...@gmail.com wrote:
> In the step-20, the mixed finite element method is used to solve a 2D
> Laplace problem, where RT and DG elements are used for the velocity and
> the pressure respectively. How does it work on spaces of Pn and Pn-1
> elements?
>
> Since the velocity is a vector, the main problem becomes representing
> the two components of the velocity both by Pn elements.

Jie Liu -- I don't think I quite understand your question. Are you
asking whether the program would work if one used Pn/Pn-1 elements
instead of RT/DG elements?

The short answer is that that's not going to work because Pn/Pn-1 are
triangular elements. But you could use Qn/Qn-1 elements. I think that's
what step-43 does, for example.

Best
W.


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

[jieli...@gmail.com]

Thank you very much for your reply. Of course, I meant Qn/Qn-1 not the FE version on triaqngles.

I have now changed the constructor as follows (following step-43) (lines 290-299):

template <int dim>
  MixedLaplaceProblem<dim>::MixedLaplaceProblem (const unsigned int degree)
    :
    degree (degree),
    //fe (FE_RaviartThomas<dim>(degree), 1,
    //    FE_DGQ<dim>(degree-2), 1),
    fe (FE_Q<dim>(degree), dim,
        FE_Q<dim>(degree-1), 1),
    dof_handler (triangulation)
  {}

The code compiles but when I run it (dim=2, degree=2 and same with degree=3), it produces the following error:

Number of active cells: 64
Total number of cells: 85
Number of degrees of freedom: 1539 (625+289)

--------------------------------------------------------
An error occurred in line <65> of file </mnt/dutita4/dv/jliu/Localdisk/dealii-8.5.0/source/dofs/dof_tools_sparsity.cc> in function
    void dealii::DoFTools::make_sparsity_pattern(const DoFHandlerType&, SparsityPatternType&, const dealii::ConstraintMatrix&, bool, dealii::types::subdomain_id) [with DoFHandlerType = dealii::DoFHandler<2>; SparsityPatternType = dealii::BlockDynamicSparsityPattern; dealii::types::subdomain_id = unsigned int]
The violated condition was:
    sparsity.n_rows() == n_dofs
Additional information:
    Dimension 914 not equal to 1539.

Stacktrace:
-----------
#0  /mnt/dutita4/dv/jliu/Localdisk/dealii-8.5.0/lib/libdeal_II.g.so.8.5.0: void dealii::DoFTools::make_sparsity_pattern<dealii::DoFHandler<2, 2>, dealii::BlockDynamicSparsityPattern>(dealii::DoFHandler<2, 2> const&, dealii::BlockDynamicSparsityPattern&, dealii::ConstraintMatrix const&, bool, unsigned int)
#1  ./step-20: Step20::MixedLaplaceProblem<2>::make_grid_and_dofs()
#2  ./step-20: Step20::MixedLaplaceProblem<2>::run()
#3  ./step-20: main
--------------------------------------------------------

I guess that the problem is here:

   

BlockDynamicSparsityPattern dsp(2, 2);
dsp.block(0, 0).reinit (n_u, n_u);
dsp.block(1, 0).reinit (n_p, n_u);
dsp.block(0, 1).reinit (n_u, n_p);
dsp.block(1, 1).reinit (n_p, n_p);
dsp.collect_sizes ();
DoFTools::make_sparsity_pattern (dof_handler, dsp);

The system matrix needs to have 3x3 blocks since the block for (n_u, n_u) no longer corresponds to a vector-valued FE space (RT) but to Qn x Qn. Hence, I have to manually generate blocks of the form (n_u_x,n_u_x), (n_u_x, n_u_y), etc. This would mean that I have to rewrite many parts of the Schur complement solver, which hard codes the 2x2 block structure of the matrix.

Now my question. Does deal.ii allow to consider Qn x Qn as an "inner" vector-valued FE space and inject this as replacement for RTn, thereby ending up in blocks (n_u, n_u), (n_u, n_p), and (n_p, n_u) and use Qn-1 as replacement for DG ending up in blocks (n_u, n_p), (n_p, n_u), and (n_p, n_p)? If this is possible, it would allow me to keep the Schur complent implementation largely unchanged, isn't this right?

Many thanks for your help in advance.

Best,

Jie Liu

[Timo Heister]

> Now my question. Does deal.ii allow to consider Qn x Qn as an "inner"
> vector-valued FE space and inject this as replacement for RTn, thereby
> ending up in blocks (n_u, n_u), (n_u, n_p), and (n_p, n_u) and use Qn-1 as
> replacement for DG ending up in blocks (n_u, n_p), (n_p, n_u), and (n_p,
> n_p)? If this is possible, it would allow me to keep the Schur complent
> implementation largely unchanged, isn't this right?

Yes. The choice on how you "block" things is independent of your FE
spaces. See step-22 for an example.

[Wolfgang Bangerth]

On 09/29/2017 06:43 AM, jieli...@gmail.com wrote:
>
> The system matrix needs to have 3x3 blocks since the block for (n_u, n_u) no
> longer corresponds to a vector-valued FE space (RT) but to Qn x Qn. Hence, I
> have to manually generate blocks of the form (n_u_x,n_u_x), (n_u_x, n_u_y),
> etc. This would mean that I have to rewrite many parts of the Schur complement
> solver, which hard codes the 2x2 block structure of the matrix.

The issue is how you compute n_u. As Timo already says, how you block things
is your choice. The problem is just that in step-20, n_u needs to be computed
differently from step-22 because in the former case, the finite element used
for the velocity is not primitive whereas in step-22 it is.