Behavior of make_hanging_node_constraints with FESystem consisting of continuous and discontinuous basis functions

[Prashant Mital]

Hello all,

I am a graduate student and I have been using deal.ii to demonstrate the use of a novel finite element method. This method has in its basis, FE_Q(1) and FE_DGQ(0) elements and I am using an FESystem object for representation. I also have local mesh refinement in parts of the domain which is why I was using make_hanging_nodes_constraints function. Once I create these hanging nodes, and try to write to the system matrix I often end up trying to access a constrained DOF and get a runtime error as a result. I was curious to know what the behaviour of the make_hanging_nodes_constraints() function will be on a system like this - the behaviour with nodal basis functions that have DOFs at vertices (FE_Q(1)) is straightforward to understand, but I cannot understand whether the function would mess with the piecewise constant DOFs (FE_DGQ(0)) and if so, how?! Furthermore, will there be any effect of the coupling DOFs between these two basis functions?

Please note that I am using make_flux_sparsity_pattern to create the system matrix and all couplings (in-cell and fluxes) are enabled.

I would appreciate any help that I can get! I am sure more experienced users have run  into this situation before.

Prashant

[Timo Heister]

Hey,

FE_DGQ has no DoFs on faces, so make_hanging_nodes_constraints will
not generate any constraints for this. So, what you are trying should
"just work". I have used FE_Q(2)^d - FE_DGQ(0) for Stokes with hanging
nodes successfully in the past. Are you using the constraint matrix
for creating the sparsity pattern and during assembly?

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

[Prashant Mital]

Dear Timo,

Thank you for your response. I will try to explain in more detail what I am trying to do.
This is what I do with the constraint matrix to create the sparsity pattern:

  constraints.clear();
  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
  constraints.close();
  BlockCompressedSparsityPattern csp(2,2);
  csp.block(0,0).reinit(n_p_cg, n_p_cg);
  csp.block(0,1).reinit(n_p_cg, n_p_dg);
  csp.block(1,1).reinit(n_p_dg, n_p_dg);
  csp.block(1,0).reinit(n_p_dg, n_p_cg);
  csp.collect_sizes();
  DoFTools::make_flux_sparsity_pattern (dof_handler, csp, constraints, true);
  sparsity_pattern.copy_from(csp);
  system_matrix.reinit (sparsity_pattern);

The method I am implementing is basically a DG method and needs the assembly of face terms which is done as per Step-30. Therefore the system matrix is being populated as follows:

for (unsigned int i=0; i<dofs_per_cell; ++i)
                          for (unsigned int j=0; j<dofs_per_cell; ++j)
                          {  system_matrix.add (dofs[i], dofs_neighbor[j],
                                      ue_vi_matrix(i,j));
                              system_matrix.add (dofs_neighbor[i], dofs[j],
                                      ui_ve_matrix(i,j));
                              system_matrix.add (dofs_neighbor[i], dofs_neighbor[j],
                                      ue_ve_matrix(i,j));
                          }

I was trying to find a way to do this using the constraint matrix class and came up with the following:

                              constraints.distribute_local_to_global(ue_vi_matrix,
                                      dofs, dofs_neighbor, system_matrix);
                              constraints.distribute_local_to_global(ui_ve_matrix,
                                      dofs_neighbor, dofs, system_matrix);
                              constraints.distribute_local_to_global(ue_ve_matrix,
                                      dofs_neighbor, dofs_neighbor, system_matrix);

I am not sure if the behaviour in this case is going to be the same as what I intended. I am hoping you can shed some light on this? Is this the right function to use?
The behaviour of the program is as follows:

1. With hanging_node_constraints and using loops to populate system matrix: Program fails after an attempt is made to access an invalid index
2. Without hanging_node_constraints and using loops to populate system matrix: Program executes and system gets solved. Convergence analysis reveals incorrect solutions
3. With hanging_node_constraints and using distribute_local_to_global to populate system matrix: UMFPACK exits with status code 1, i.e. system_matrix turns out to be singular.

Does this seem correct? What could be causing the matrix to become singular in case 3?

[Timo Heister]

> I was trying to find a way to do this using the constraint matrix class and
> came up with the following:

Option A:
assemble manually, use constraints.condense(matrix, vector);

Option B:
assemble with constraints.distribute_local_to_global and do not condense.

> Does this seem correct? What could be causing the matrix to become singular
> in case 3?

Note sure, but I would recommend option B. Does it work if you don't
have hanging nodes?

[Prashant Mital]

Okay so using option B gives the correct results when I run the code with make_hanging_node_constraints but the mesh is matching and has no hanging nodes/local refinement. I then ran the case with one cell refined and UMFPACK starts returning the same error code pertaining to my matrix being singular. So I suppose one possibility is that I am handling non-matching interfaces improperly and I am going to go back and check that. What are the other possible sources of error in this situation as far the constraints are concerned  (in conjunction with the FE_Q(1) and FE_DGQ(0) space). Do constraints affect what happens to the off diagonal terms in the system matrix (e.g. the cross-terms between the two types of basis functions)?

[Prashant Mital]

So the method I am using has a weak form that loops over cells and interfaces (as opposed to a weak form that is integrated over the entire domain). As a result, I think that the system does not need any special treatment of hanging nodes using the Constraint Matrix class.

Timo, you mentioned that you used this type of a basis for Stokes flow - would you mind sharing what your weak form looked like? Was it akin to a continuous galerkin setup with integrals set up over the entire solution domain, or local integrals as is seen in discontinuous galerkin methods? If it was the latter, I have a feeling the constraints might not be necessary.

I will try to redo the convergence analysis and update here in case I can get the method to work.

Thanks for the help!
Prashant

[Timo Heister]

> What are the other possible sources of error in this situation as far the constraints are concerned (in conjunction with the FE_Q(1) and FE_DGQ(0) space).

Are you using the constraint matrix whenever you enter stuff into
matrices/vectors? Are you distributing constraints after solving and
before you use the solutions for postprocessing/output/following
solution steps?

> Do constraints affect what happens to the off diagonal terms in the system matrix (e.g. the cross-terms between the two types of basis functions)?

Not sure how that could be a problem.

> Timo, you mentioned that you used this type of a basis for Stokes flow -
> would you mind sharing what your weak form looked like? Was it akin to a
> continuous galerkin setup with integrals set up over the entire solution
> domain, or local integrals as is seen in discontinuous galerkin methods? If
> it was the latter, I have a feeling the constraints might not be necessary.

Standard continuous Galerkin of Stokes.

> I will try to redo the convergence analysis and update here in case I can
> get the method to work.

Do you have face integrals? How are you handling them when you have a
face that is refined from one side?

[Prashant Mital]

Dear Timo,

Lack of theoretical background/literature in what I am trying to achieve seems to be hampering some progress as far as checking convergence rates goes. However, I do think that constraints are not required to be imposed in my method as the hanging nodes dont need any special treatment. This is because my test space has a component which has only local support and so all my integrals are locally defined. Not sure yet how that works mathematically which is what needs more scrutiny. 

I do have face integrals. I only operate on a face if I approach it from the coarser side. Once I am there, I reinit FE(Sub)FaceValues objects using both neighbouring cells common to a subface and then I perform my computations. I use 4 matrices, ui_vi, ui_ve_ue_vi and ue_ve, one a subface has been operated on, I used the constraints to add the info in the global matrix and move onto the next subface. As I said, the implementation is algorithmically identical to step-31.