Error with Simplex and adaptive mesh with FESubfaceValues

[Sanghyun Lee]

Hello!

Happy New Year. 

I am trying to implement 

- adaptive mesh refinement (with hanging nodes) for an elasticity equation 

- with DG methods 

- by utilizing Simplex. 

I am willing to modify library files if necessary. 

Please if there is anybody who has looked through above components, please share your thoughts! :) 

For DG methods, I am using the classical face loop, so we have three different cases. 

When the face has children (neighbor is finer)-   we need to use "FESubfaceValues".. 

---

    MappingFE<dim>     mapping;    mapping(FE_SimplexP<dim>(1)

    FESystem<dim> fe; fe(FE_SimplexDGP<dim>(1), dim)
    QGaussSimplex<dim> face_quadrature_formula(fe.degree+1); 

    

     FESubfaceValues<dim> fe_subface_values (mapping,fe, face_quadrature_formula,
       update_values    | update_normal_vectors |
       update_gradients  |
        update_quadrature_points  | update_JxW_values);

---

then, we initialize with cell, face and subface numbers - 

     fe_subface_values.reinit(cell,face_no,subface_no);

However, with the simplex, I have the following error with the above line ! 

An error occurred in line <1597> of file </dealii/source/base/qprojector.cc> in function

static QProjector<2>::DataSetDescriptor dealii::QProjector<2>::DataSetDescriptor::subface(const dealii::ReferenceCell &, const unsigned int, const unsigned int, const unsigned char, const unsigned int, const internal::SubfaceCase<2>) [dim = 2]

The violated condition was: 

reference_cell == ReferenceCells::Quadrilateral

My question is 

i) does anybody have any experience with this case? when you need subface values for DG + Simplex? 

ii) Should I use meshworker or something more recent? 

iii) I looked into qprojector.cc, and QProjector<2>::DataSetDescriptor::subface.

This function only returns 

return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
          n_quadrature_points);

If I manually fix this part to fit simplex, and reinstall deal.ii, do you think it will fix the problem? 

Thanks a lot! 

Sanghyun

[Peter Munch]

Hi Sanghyun,

there is an open PR that should (hopefully) fix the issue: https://github.com/dealii/dealii/pull/17908. It is not merged yet. But you could try it out and gives us feedback!

Best,
Peter

[Sanghyun Lee]

Dear Peter

Thanks for the reply and pointing the direction that I was exactly looking for. 

I downloaded and installed with the modified files on the PR, and checked with my code. The initial concern is now removed.  

However, now I found one more discrepancy for the simplex case in the following code line: 

Assert (neighbor->neighbor_child_on_subface(neighbor_face_subface.first, neighbor_face_subface.second)
== cell, ExcInternalError());

Please see the figure:

1. Let's assume that we are in the [cell 1] and  we are on the blue face between [cell 1] and [cell 0]

2. Then, the neighbor of [cell 1] (across blue face) is [cell 0], which is coarser.

Now, cell = cell 1, neighbor = cell 0.

3. In this case, the above code has an issue. 

Why? 

The output of 

"neighbor->neighbor_child_on_subface(neighbor_face_subface.first, neighbor_face_subface.second)"

will give me [cell 2]. (note that neighbor_face_subface.first = 1,  neighbor_face_subface.second = 0, which makes sense)

It should give me [cell 0]! 

But, it is interesting that in the Simplex case, [cell 0] is thinking [cell 2] is also the neighbor->child !? (I think I know this could be a case for Simplex...)  

4. Still, note that it is interesting the output of "neighbor->face(neighbor_face_subface.first)->n_children()" 

will give me still 2 (which is good). 

I am trying to figure this out, but if you have any idea, please let us know. 

Sanghyun

[Sanghyun Lee]

I fixed the issue; but before I discuss that, 

=========

first, sorry there was a typo in my above email. The correct discussion is

... 

Assert (neighbor->neighbor_child_on_subface(neighbor_face_subface.first, neighbor_face_subface.second)
== cell, ExcInternalError());

Please see the figure:

1. Let's assume that we are in the [cell 1] and  we are on the blue face between [cell 1] and [cell 0]

2. Then, the neighbor of [cell 1] (across blue face) is [cell 0], which is coarser.

Now, cell = cell 1, neighbor = cell 0.

3. In this case, the above code has an issue. 

The output of 

"neighbor->neighbor_child_on_subface(neighbor_face_subface.first, neighbor_face_subface.second)"

will give me [cell 3]. (note that neighbor_face_subface.first = 1,  neighbor_face_subface.second = 0, which makes sense)

It should give me [cell 1]! 

=========

Now, here is how I fixed the issue. 

- In source/grid/tria_accessor.cc

- in the function "neighbor_child_on_subface"

I changed 

the following line 

const unsigned int neighbor_subface = (!(this->line_orientation(face)) != !(neighbor_cell->line_orientation(neighbor_face))) ? (1 - subface) : subface;

to

if(neighbor_cell->face_orientation(neighbor_face) == neighbor_cell->face_orientation(face)) 

{ 

   if(neighbor_cell->face_orientation(neighbor_face) == 0) 

           neighbor_subface = 1 - subface; 

  else 

            neighbor_subface = subface; 

} 

else

{ 

    if(neighbor_cell->face_orientation(neighbor_face) == 0) 

           neighbor_subface = 1 - subface; 

    else 

           neighbor_subface = subface; 

}

Why? this is due to the orientation and numbering of the faces for the Simplex (I think). 

To be honest, I am not sure this is the right way to do it or efficient or general way to do it or not. At least for now, I do have expected 'convergence rates' for 'my' problem. 

If anyone thinks this is not correct or has a better idea, please share. 

Thanks

Sanghyun

[Wolfgang Bangerth]

On 1/8/25 22:00, Sanghyun Lee wrote:
>
> *Now, here is how I fixed the issue. *
> [...]

Would you mind submitting a pull request with your changes? That's the easiest
way to discuss the correctness of what you have.

Best
W.

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