Method to find an "eggshell" subdomain

[Lukas Korous]

Hello,

we want to implement the method described in the attached paper in deal.ii. Basically for an area with a marker "m" and for a given parameter p (integer) we need to find all elements that are <=p elements far from the border of area "m" without those in area "m" - this is meant by the word "eggshell".

Plus we need to:

1) identify the two boundaries of such a set of elements (in case "m" is a circle, the eggshell will be an annulus)
2) set an FE space on this mesh
3) set Dirichlet conditions to the on the two boundaries from 1)
4) calculate a simple equation, with solution in the FE space from 2) respecting BCs from 3)

But this will probably be easy, I would like to ask for a kick in the right direction how to find the right elements.

Thank you,
Lukas

[Wolfgang Bangerth]

Lukas,

> we want to implement the method described in the attached paper in deal.ii.
> Basically for an area with a marker "m" and for a given parameter p (integer)
> we need to find all elements that are <=p elements far from the border of area
> "m" without those in area "m" - this is meant by the word "eggshell".

Let's say you marked the cells in your area by setting its user index to 1,
the remaining cells having user index 0. Then you'd want something like this:

for(unsigned int level=1; level<=p; ++level)
for(cell=...)
if (cell->user_index() == level) // cell was tagged in the last round
for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if(cell->at_boundary(f) == false)
if(cell->neighbor(f)->user_index()==0) // we found a neighbor that
// was not previously tagged
cell->neighbor(f)->set_user_index(level+1);

This way, you will tag all cells that are face neighbors up to 'p' levels away
with the level number (from 2 to p+1). It's a bit more complicated if you have
an adaptively refined mesh, but the principle is the same. It requires a more
complicated algorithm if you want to also tag vertex neighbors or, in 3d, edge
neighbors (not just face neighbors).

> Plus we need to:
>
> 1) identify the two boundaries of such a set of elements (in case "m" is a
> circle, the eggshell will be an annulus)

These would be faces where cells of different user indices come together.

> 2) set an FE space on this mesh

You mean create a space *only* on the eggshell cells? That would be a case for
hp::DoFHandler and FENothing.

> 3) set Dirichlet conditions to the on the two boundaries from 1)

There is no function to set Dirichlet conditions in the interior of a cell
right now. But it should not be very difficult to copy-paste the existing
VectorTools::interpolate_boundary_values() function and modify it to run on
internal boundaries.

Best
W.

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

[Lukas Korous]

Hello Wolfgang,

thank you for the hints - I will try to do this.

Best,
Lukas

[Jean-Paul Pelteret]

FYI: This discussion may also be relevant to your primary problem.