Subdivided cylinder and boundary ids

[Alex Cumberworth]

I have been trying to construct a system that is small (a total length of 1e-5 m) with the subdivided cylinder grid generator. However, I have been having trouble with the boundary indicators. I wanted to set my own boundaries, and because there is no colorize option (it defaults to setting the left end face of the cylinder to 1, and the right end to 2), I tried to first set all boundary ids to 0. However, this leads to the following exception:

terminate called after throwing an instance of 'dealii::TriaAccessorExceptions::ExcNoPeriodicNeighbor'
  what():  
--------------------------------------------------------
An error occurred in line <2597> of file </home/ipausers/cumberworth/src/dealii-9.3.0/source/grid/tria_accessor.cc> in function
    dealii::TriaIterator<dealii::CellAccessor<dim, spacedim> > dealii::CellAccessor<dim, spacedim>::neighbor_or_periodic_neighbor(unsigned int) const [with int dim = 3; int spacedim = 3]
The violated condition was:
    false
Additional information:
    (none)

Since in my test system, I was setting the boundary ids to be the same as the default, I tried to rely on this instead. However, I was having strange behaviour where it seemed that the entire system was having the boundary conditions applied to it. I scrolled through grid_generator.cc and found that subdivided_cylinder has a hardcoded epsilon of 1e-5 for determining whether a face is at the boundary. Of course, with a system of my size, all faces will end up being considered to be at the boundary.

Because it not only sets the boundary ids, but the manifold ids, it seems I will have to modify the subdivided_cylinder to make the epsilon value much smaller or somehow relative to the system size, which is simple enough.

It would probably also be good to have the colorize option. However, I am not sure what to make of the exception when I try to set all the boundary ids to 0.

Best,

Alex

[Marc Fehling]

Hi Alex,

> I scrolled through grid_generator.cc and found that subdivided_cylinder has a hardcoded epsilon of 1e-5 for determining whether a face is at the boundary. Of course, with a system of my size, all faces will end up being considered to be at the boundary.

>

> Because it not only sets the boundary ids, but the manifold ids, it seems I will have to modify the subdivided_cylinder to make the epsilon value much smaller or somehow relative to the system size, which is simple enough.

That sounds like a good idea to me. Would you willing to write a patch and open a pull request for the library? I'm sure others will benefit of your changes as well!

There are several options you can consider.

1. You could provide epsilon as a parameter to the function and default it to `1e-5`.
   
2. Instead of an absolute threshold, you can introduce a relative one, for example, `epsilon * cell->diameter()` or something similar.
3. You can also try to use the machine precision std::numeric_limits<T>::epsilon as a substitute for `1e-5`.
   
4. Why do we need these epsilons in the first place? Maybe it will work if you just replace the greater/lesser checks with their equal variants `>=` and `<=`.
   

Best,
Marc

[Marc Fehling]

By the way, this function has just been introduced in April with #12003. Suggestions for improvement on this new feature are welcome!

Best,
Marc

[Wolfgang Bangerth]

On 6/25/21 12:37 PM, Marc Fehling wrote:
> # Instead of an absolute threshold, you can introduce a relative one, for

> example, `epsilon * cell->diameter()` or something similar.

FWIW, in most places, this is what we do.
Best
W.


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

[blais...@gmail.com]

I most likely copied this mistake from the regular Cylinder and did not think of it when I coded the subdivided cylinder.
Alex, if you feel comfortable, you could go ahead and make a PR following Wolfgang's suggestion. Otherwise, I'll gladly look into it this week :)

Best

BRuno

[Alex Cumberworth]

Hi Bruno,

Sure, I can try making a PR this week.

Best,

Alex

[blais...@gmail.com]

This issue is now fixed by a PR merged yesterday.

Best

Bruno