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 *neighbo*r 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
On Wednesday, January 8, 2025 at 8:00:39 PM UTC-5 Sanghyun Lee wrote:
> 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 *neighbo*r 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
>
> [image: ParaView_5_11_1-2.jpg]
>
>
> On Wednesday, January 8, 2025 at 3:23:32 PM UTC-5 peterr...@gmail.com
> wrote:
>
>> 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
>>
>> On Wednesday, 8 January 2025 at 21:10:08 UTC+1 shl....@gmail.com wrote:
>>
>>> 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
>>>
>>>
>>>
--
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+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/dealii/e8663f5f-a2e7-4307-bcae-2beaa47ea1den%40googlegroups.com.
