Error in step-51 when using 3-dim with adaptivity

[Jau-Uei Chen]

Dear dealii community,

I am currently working on a HDG solver and use step-51 as a reference. I found that the example does not work for a 3d setting along with adaptivity. It seems that an error is from "VectorTools::project_boundary_values", the error message I got is:

An error occurred in line <742> of file </org/groups/ceo/chenju/candi/install/tmp/unpack/deal.II-v9.2.0/include/deal.II/numerics/vector_tools_boundary.templates.h> in function                                                                                            
    void dealii::VectorTools::internal::do_project_boundary_values(const M_or_MC<dim, spacedim>&, const DoFHandlerType<dim, spacedim>&, const std::map<unsigned int, const dealii::Function<spacedim, number>*>&, const Q_or_QC<(dim - 1)>&,std::map<unsigned int, number>&, std::vector<unsigned int>) [with int dim = 3; int spacedim = 3; DoFHandlerType = dealii::DoFHandler; M_or_MC = dealii::Mapping; Q_or_QC = dealii::Quadrature; number = double]                                                                                            
The violated condition was:                                                                                                          
    level == cell->level()                                                                                                            
Additional information:                                                                                                              
    The mesh you use in projecting boundary values has hanging nodes at the boundary. This would require dealing with hanging node constraints when solving the linear system on the boundary, but this is not currently implemented. 

The version of dealii I use is 9.2.0. I am wondering if there is any way to resolve this error. Besides this error, I also have a hard time investigating the vertice initiated by FE_FaceQ and would also like to seek some suggestions. Here is my code modified from step-51:

template <int dim>
  void HDG<dim>::assemble_system_one_cell(
    const typename DoFHandler<dim>::active_cell_iterator &cell,
    ScratchData &                                         scratch,
    PerTaskData &                                         task_data)
  {

       for( unsigned int face = 0; face < GeometryInfo< dim >::faces_per_cell; ++face ){

            if(cell->face(face)->at_boundary() && (cell->face(face)->boundary_id() == 0)){

                 std::cout << "node points" << std::endl;
                 for (const auto v: cell->face(face)->vertex_indices()){

                      const Point<dim> node_point = cell->face(face)->vertex(v);
                       std::cout << node_point << std::endl;
                  }

            }

       }

  }

But got an error saying:

‘const class dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false>’ has no member named ‘vertex_indices’; did you mean ‘vertex_index’?  

Any suggestion or comment is greatly appreciated.

Best Regards,

Jau-Uei Chen

[Wolfgang Bangerth]

> I am currently working on a HDG solver and use step-51 as a reference. I found
> that the example does not work for a 3d setting along with adaptivity. It
> seems that an error is from "VectorTools::project_boundary_values", the error
> message I got is:
>
> An error occurred in line <742> of file
> </org/groups/ceo/chenju/candi/install/tmp/unpack/deal.II-v9.2.0/include/deal.II/numerics/vector_tools_boundary.templates.h>
> in function
>     void dealii::VectorTools::internal::do_project_boundary_values(const
> M_or_MC<dim, spacedim>&, const DoFHandlerType<dim, spacedim>&, const
> std::map<unsigned int, const dealii::Function<spacedim, number>*>&, const
> Q_or_QC<(dim - 1)>&,std::map<unsigned int, number>&, std::vector<unsigned
> int>) [with int dim = 3; int spacedim = 3; DoFHandlerType =
> dealii::DoFHandler; M_or_MC = dealii::Mapping; Q_or_QC = dealii::Quadrature;
> number = double]
> The violated condition was:
>     level == cell->level()
> Additional information:
>     The mesh you use in projecting boundary values has hanging nodes at the
> boundary. This would require dealing with hanging node constraints when
> solving the linear system on the boundary, but this is not currently implemented.
>
> The version of dealii I use is 9.2.0. I am wondering if there is any way to
> resolve this error.

Not easily. Basically, what is happening is that for *projection* of boundary
values, you need to build and solve a linear system of only those degrees of
freedom that are located on the boundary. The way this is currently done is
unable to deal with hanging node constraints.

There are approaches to deal with this, but they are (as far as I know) not
currently implemented and so even newer versions of deal.II will not support
this. There are two options for you:
* You can use *interpolation* (rather than projection) of boundary values
* You can implement the missing functionality. We'd be happy to guide you on
how this can be done, but it will involve a certain amount of work.

> Besides this error, I also have a hard time investigating
> the vertice initiated by FE_FaceQ and would also like to seek some
> suggestions. Here is my code modified from step-51:
>
> template <int dim>
>   void HDG<dim>::assemble_system_one_cell(
>     const typename DoFHandler<dim>::active_cell_iterator &cell,
>     ScratchData &                                         scratch,
>     PerTaskData &                                         task_data)
>   {
>        for( unsigned int face = 0; face < GeometryInfo< dim
> >::faces_per_cell; ++face ){
>             if(cell->face(face)->at_boundary() &&
> (cell->face(face)->boundary_id() == 0)){
>                  std::cout << "node points" << std::endl;
>                  for (const auto v: cell->face(face)->vertex_indices()){

This uses a function that was introduced only in deal.II 9.3, I believe. But
you can equivalently implement this line using
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)


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

[Jau-Uei Chen]

Dear Dr. Wolfgang Bangerth,

Thanks for your response. I am quite interested in how to implement projection. 

In fact, I have tried both methods but did not get the correct result using projection:

• Interpolate the boundary data: It works quite well for adaptivity with both 2 and 3 dimensions. Though expecting extra error will be introduced, I observed that it is actually comparable to the case with projected boundary data (I do the experiment through global refinement).
• Project the boundary data:  I tried to project the boundary data by commenting out the project_boundary_values and adding attached code in step-51:

    if (task_data.trace_reconstruct == false)
      {
        ...
        cell->get_dof_indices(task_data.dof_indices);

       my code 

     }else{

       ...

     }

where face_mass is the size of fe.dofs_per_cell x fe.dofs_per_cell matrix, and both dirichlet_rhs and proj_dirichlet_rhs are fe.dofs_per_cell vector. I just solve the linear system on the boundary faces, assign the projected data to the cell_vector, and zero out the off-diagonal terms. In this procedure, I assume constraints can handle hanging nodes properly. Whatever values I assign to the corresponding row will be replaced by some algebraic function that will be used to maintain the continuity via data from end support points.

The code can run but the linear solver for solving trace unknowns requires a significant iteration number (will explode if using global refinement option). In addition, the L2 error norm is too high compared to the original setting. Hence, I may do something wrong. Could you give me some suggestions or comments on it?

Best Regards,

Jau-Uei Chen

Wolfgang Bangerth 在 2022年3月20日 星期日下午9:49:49 [UTC-5] 的信中寫道：

[Wolfgang Bangerth]

> In fact, I have tried both methods but did not get the correct result using
> projection:
>

> * Interpolate the boundary data: It works quite well for adaptivity with

> both 2 and 3 dimensions. Though expecting extra error will be introduced,
> I observed that it is actually comparable to the case with projected
> boundary data (I do the experiment through global refinement).

Yes, in practice there is little difference between the accuracy obtained from
using boundary value interpolation and projection. Because the former is
computationally cheaper, that's what everyone uses.

If it works for you, why not go with it?


> * Project the boundary data:  I tried to project the boundary data by

> commenting out the project_boundary_values and adding attached code in
> step-51:
>
>     if (task_data.trace_reconstruct == false)
>       {
>         ...
>         cell->get_dof_indices(task_data.dof_indices);
> my code
>      }else{
>        ...
>      }
>
> where face_mass is the size of fe.dofs_per_cell x fe.dofs_per_cell matrix, and
> both dirichlet_rhs and proj_dirichlet_rhs are fe.dofs_per_cell vector. I just
> solve the linear system on the boundary faces, assign the projected data to
> the cell_vector, and zero out the off-diagonal terms. In this procedure, I
> assume constraints can handle hanging nodes properly. Whatever values I assign
> to the corresponding row will be replaced by some algebraic function that will
> be used to maintain the continuity via data from end support points.
>
> The code can run but the linear solver for solving trace unknowns requires a
> significant iteration number (will explode if using global refinement option).
> In addition, the L2 error norm is too high compared to the original setting.
> Hence, I may do something wrong. Could you give me some suggestions or
> comments on it?

You'll have to debug it, I'm afraid :-) Make sure it is correct first, for
example by imposing boundary conditions equal to 1.0 and ensuring that what
you get on each cell/face is correct. You can test this on a cube with only
one cell first, then refine once, etc. In your code, I don't understand why
you need to multiply
scratch.proj_dirichlet_rhs(ii) * task_data.cell_matrix(ii,ii);
when assigning to cell_vector. Is that because you want to set the
corresponding entries in the rhs as boundary values? If so, is the cell_matrix
already correctly computed?

If I understand your approach correctly, then you are doing a local projection
on each face, one face after the other. You will then overwrite the already
computed values of DoFs on a face that are shared with a previously visited
face. That's ok, but it's not the same as a global projection, of course. I
suspect that this will not be more accurate than interpolation of boundary values.

[Jau-Uei Chen]

Thanks for Dr. Wolfgang Bangerth's detailed information and very great tip on debugging though I am settled with interpolation now. However, I would like to seek some suggestions about implementation on FE_FaceQ defined on skeleton that conforming to the small elements.

scratch.proj_dirichlet_rhs(ii) * task_data.cell_matrix(ii,ii);
when assigning to cell_vector. Is that because you want to set the
corresponding entries in the rhs as boundary values? If so, is the cell_matrix
already correctly computed?

Yes, I want to set the corresponding entries in the rhs as boundary values. Since the implementation is placed in the end of  HDG::assemble_system_one_cell, I believe that the cell_matrix is already computed. 

If I understand your approach correctly, then you are doing a local projection
on each face, one face after the other. You will then overwrite the already
computed values of DoFs on a face that are shared with a previously visited
face. That's ok, but it's not the same as a global projection, of course. I
suspect that this will not be more accurate than interpolation of boundary values.

Thanks for the pointing it out. When I implement it, I thought that at boundary we have faces conforming to the fine elements since there is no adjacent element. However, I just realized that it is not the case for dealii when I print out the information about DoFs on a simple nonconforming mesh.

As pointed out at beginning,  I am also wondering if it is possible to implement FE_FaceQ in a way that it is defined on skeleton that conforming to the small elements. That is, no constraint is required.  According to the suggestions given in the posts "Fe_faceq on adapted grids" and "Error when make_hanging_node_constraints: FE_TraceQ element", it seems that it is possible though no one actually try to do it. 

Best Regards,

Jau-Uei Chen

 

[Wolfgang Bangerth]

On 3/31/22 09:42, Jau-Uei Chen wrote:
>
> As pointed out at beginning,  I am also wondering if it is possible to
> implement FE_FaceQ in a way that it is defined on skeleton that conforming to
> the small elements. That is, no constraint is required.

I don't quite understand what you are trying to do. Can you explain in more
detail?

[Jau-Uei Chen]

Sorry for causing the confusion. What I am looking for is to implement HDG with the discontinuous finite element defined in the splitted mortars (blue):

 0 - - - - - -  0 0  0 - - 0 0 - - 0 

 |                |   |   |      |  |       |

 |                |  0  0 - - 0 0 - - 0       

 |                |  0  0 - - 0 0 - - 0  

 |                |   |   |      |  |       |    

 0 - - - - - -  0 0  0 - - 0 0 - - 0

instead of the full-sided (red) one which is currently used in dealii where we put the constraint and hence get a full-sided one. 

 0 - - - - - -  0 0  0 - - 0 0 - - 0 

 |                |   |   |      |  |       |

 |                |   |  0 - - 0 0 - - 0    

 |                |  0   

 |                |   |  0 - - 0 0 - - 0  

 |                |   |   |      |  |       |    

 0 - - - - - -  0 0  0 - - 0 0 - - 0

Is it possible to do so in dealii?

Best Regards,

Jau-Uei Chen

Wolfgang Bangerth 在 2022年3月31日 星期四下午3:00:30 [UTC-5] 的信中寫道：

[Wolfgang Bangerth]

On 3/31/22 14:36, Jau-Uei Chen wrote:
> *** Caution: EXTERNAL Sender ***

>
> Sorry for causing the confusion. What I am looking for is to implement
> HDG with the discontinuous finite element defined in the splitted
> mortars (blue):
>
>  0 - - - - - -  0 0  0 - - 0 0 - - 0
>  |                | |   |      |  |       |
>  |                | 0  0 - - 0 0 - - 0
>  |                | 0 0 - - 0 0 - - 0
>  |                | |   |      |  |       |
>  0 - - - - - -  0 0  0 - - 0 0 - - 0
>
> instead of the full-sided (red) one which is currently used in dealii
> where we put the constraint and hence get a full-sided one.
>
>  0 - - - - - -  0 0  0 - - 0 0 - - 0
>  |                | |   |      |  |       |
>  |                | |  0 - - 0 0 - - 0
>  |                | 0
>  |                | |  0 - - 0 0 - - 0
>  |                |  |   |      |  |       |
>  0 - - - - - -  0 0  0 - - 0 0 - - 0
>
> Is it possible to do so in dealii?

(A comment first: I think in the second picture, you wanted to use only
two red DoFs, not three, on the common edge, right?)

I see. That may be possible, but it is difficult.

The way things are implemented is that there is a finite element space
associated with every edge in the mesh (in 2d). This means that there is
a finite element space associated with the coarse edge of the left cell,
and one with each of the small edges of the cells on the right, and then
constraints are used to make sure that the two are the same.

What you now want is to associate a finite element space with the coarse
face *that depends on the fact that that face is refined*. In other
words, you would want to use a different finite element space for the
four edges of the cell on the left: One space for the left, top, and
bottom face, and a different one for the face on the right that has
neighboring cells. This is not how the internal data structures are set
up currently.

If you *had* to implement something like this, then the way to do it
would probably be to implement the space that corresponds to your top
example and use it on *all* edges. Then you would create hanging node
constraints as usual, and after that you would loop over the mesh and
add constraints yourself for all edges that are not further refined;
these constraints would bring the 2+2 DoFs back down to the two you want
on these edges.

This approach is not convenient, because it's not already implemented in
deal.II, but it probably wouldn't be terribly hard to implement. You'd
have to re-implement something like FEFaceQ, though, because you now
want the space to not be defined on a face, but on each of its children
(whether the edge is refined or not doesn't actually matter here). But
it can be done, and you might be interested in seeing something like
this in practice in this paper:

https://www.math.colostate.edu/~bangerth/publications/2016-nonconforming.pdf

[Jau-Uei Chen]

Thanks for the detailed explanation. It is quite helpful. However, I am wondering if you could elaborate on implementation side a little bit more. 

(A comment first: I think in the second picture, you wanted to use only
two red DoFs, not three, on the common edge, right?)

 

Yes, that is right and I believe that it is how dealii originally implement FE_FaceQ after enforcing the constraints.

If you *had* to implement something like this, then the way to do it
would probably be to implement the space that corresponds to your top
example and use it on *all* edges. Then you would create hanging node
constraints as usual, and after that you would loop over the mesh and
add constraints yourself for all edges that are not further refined;
these constraints would bring the 2+2 DoFs back down to the two you want
on these edges.

This approach is not convenient, because it's not already implemented in
deal.II, but it probably wouldn't be terribly hard to implement. You'd
have to re-implement something like FEFaceQ, though, because you now
want the space to not be defined on a face, but on each of its children
(whether the edge is refined or not doesn't actually matter here). 

It is not clear to me how modifying FEFaceQ can even achieve the first step (highlighted in yellow) you describe. It seems that FEFace_Q is universally defined on *the unit face* and cannot see the structure of triangulation. Hence, the finite element space is always defined on the each face of an element when we are doing "distribute_dofs". Are you referring that we can use DSSY element (given in the paper) to construct the discontinuous function defined on the reference domain of a face to achieve the first step? 

An possible way I imagine is to re-defined FEFace_Q in the following way to complete the first step: 

(assume degree = 1, so fcn 1 and fcn 2 are  Lagrange polynomial with degree 1)

    fcn1          fcn2

*------------**------------*     <- basis function with supporting points *

o - - - - - - - - - - - - - o     <- geometric line segment in reference domain with vertices o

However, it seems not possible since the construction of a finite element space relies on inheritance of "FE_PolyFace" and such discontinuous function is not polynomial.

Any further comment or suggestion is greatly appreciated.

Best Regards,

Jau-Uei Chen

 

[Wolfgang Bangerth]

On 4/1/22 14:40, Jau-Uei Chen wrote:
>
> It is not clear to me how modifying FEFaceQ can even achieve the first
> step (highlighted in yellow) you describe. It seems that FEFace_Q is
> universally defined on *the unit face* and cannot see the structure of
> triangulation. Hence, the finite element space is always defined on the
> each face of an element when we are doing "distribute_dofs".

Correct.

> Are you
> referring that we can use DSSY element (given in the paper) to construct
> the discontinuous function defined on the reference domain of a face to
> achieve the first step?

No, what I meant is that there are other contexts in which one wants to
manually add constraints, and the situation in the linked paper is one.
In your case, you would want to add additional constraints for all faces
that are *not* in fact refined. This is because on these faces you want
to bring the space back down to linear and continuous across the entire
face.

> An possible way I imagine is to re-defined FEFace_Q in the following way
> to complete the first step:
> (assume degree = 1, so fcn 1 and fcn 2 are  Lagrange polynomial with
> degree 1)
>
>     fcn1          fcn2
> *------------**------------*     <- basis function with supporting points *
> o - - - - - - - - - - - - - o     <- geometric line segment in reference
> domain with vertices o

Yes, this is exactly how it should look.

> However, it seems not possible since the construction of a finite
> element space relies on inheritance of "FE_PolyFace" and such
> discontinuous function is not polynomial.

No, the inheritance is not necessary. FE_PolyFace is just a convenience
but not a requirement. At the end of the day, you need to implement the
interface of the FiniteElement class, and how you do this is up to you.
Different FE_* classes do this in different ways, depending on the
structure of their finite element spaces. You can definitely implement
non-polynomial finite elements -- you just need to implement all of the
required virtual functions.

(I should also mention that because your shape functions are only
piecewise polynomial, whenever you integrate terms on these faces you
need to use iterated quadrature formulas. QIterated is your friend.)