Creating an internal (curved) interface in deal.II -- a (failed) attempt and some questions

[Sean Carney]

Hi all,

I want to create a 2D domain in deal.II with an internal, curve interface $\Gamma$. A simple example of what I want to create is in the attached screenshot (credit to Bartels, Bonito & Tscherner).

In short, does anyone have any suggestions for how to do this?

The critical thing is that upon mesh refinement, deal.II knows where my interface is and adds points accordingly (of course, this issue is already discussed in Step-1). Here, however, the 1D manifold is internal, so maybe this changes things?

What I tried so far is the following: 

I have an analytic expression for my curve (i.e. it can be written in the form $(\Gamma(y), y)$), so I created two different triangulation objects using the

   GridGenerator::subdivided_hyper_rectangle

and 

   GridTools::transform

functions. Then, I tried to merge them together. However, I am receiving an error:

"An error occurred in line <1235> of file <dealii/source/grid/grid_tools.cc> in function
    void dealii::GridTools::{anonymous}::AdjacentCells<2>::push_back(const dealii::GridTools::{anonymous}::AdjacentCell&)
The violated condition was:
    adjacent_cells[1].cell_index == numbers::invalid_unsigned_int"

Of course, I know that the vertices along the common interface, in this case $\Gamma$, need to be aligned. I believe that this is indeed the case.

My code is simple enough:

void grid_10()
{
  const int N = 14;

   // create 1st triang
  Triangulation<2>          tri_1;

  std::vector<unsigned int> repetitions(2);
  repetitions[0] = N;
  repetitions[1] = N;
  GridGenerator::subdivided_hyper_rectangle(tri_1,
                                            repetitions,
                                            Point<2>(0.0, 0.0),
                                            Point<2>(2./3., 1.));

  GridTools::transform(
    [](const Point<2> &in) {
      if (in[0] < 0.5)
      {
        return in;
      }
      else
      {
        return Point<2>(2./3. - 1./6.*std::sin(numbers::PI * in[1]), in[1]);
      }
    },
    tri_1);

   // create 2nd triang
  Triangulation<2>          tri_2;
 
  repetitions[0] = N;
  repetitions[1] = N;
  GridGenerator::subdivided_hyper_rectangle(tri_2,
                                            repetitions,
                                            Point<2>(2./3., 0),
                                            Point<2>(1., 1.));
 
  GridTools::transform(
    [](const Point<2> &in) {
      if (in[0] > 0.7)
      {
        return in;
      }
      else
      {
        return Point<2>(2./3. - 1./6.*std::sin(numbers::PI * in[1]), in[1]);
      }
    },
    tri_2);

   // merge them
  Triangulation<2>          triangulation;
  GridGenerator::merge_triangulations(tri_1,tri_2,triangulation);
  print_mesh_info(triangulation, "grid.vtu");
 
}

Is there an obvious error that I am missing?

Even more critically, will this even accomplish what I want? I see in the "merge_triangulations" documentation that " In particular, manifolds attached to the two input triangulations will be lost in the result triangulation."

Any help is greatly appreciated. Thank you!

Best regards,

-Sean

[Winnifried Wollner]

Dear Sean,
I think you can find an example for what you want to achieve in the code
of a recent paper https://epubs.siam.org/doi/10.1137/23M158262X

The code is available at

https://zenodo.org/records/10459309

where you can check the file
InverseHomogenization/geometry.h

Best regards
Winni

On 05.09.24 22:13, Sean Carney wrote:
> Hi all,
>
> I want to create a 2D domain in deal.II with an internal, curve interface
> $\Gamma$. A simple example of what I want to create is in the attached
> screenshot (credit to Bartels, Bonito & Tscherner

> <https://doi.org/10.1093/imanum/drad004>).

--
--------------------------------------------------------
Prof. Dr. Winnifried Wollner
Universität Hamburg
MIN Faculty
Department of Mathematics
Bundesstr. 55; 20146 Hamburg; Germany
--------------------------------------------------------
Room : 115
Phone: +49 (0)40 42838-4079
Web : https://www.math.uni-hamburg.de/personen/wollner/
eMail: winnifrie...@uni-hamburg.de
--------------------------------------------------------

[Sean Carney]

Dear Winni,

Thank you very much for pointing me to the helpful example from the code from your recent paper.

I have learned a great deal from breaking down this function and understanding what each part is 

doing (thank you for adding the descriptive comments within).

I think I now understand reasonably well how to create a mesh with some kind of internal interface $\Gamma$

*so long as* the interface is described by one of the built-in functions within the GridGenerator namespace

(for example, hyper_ball, cylinder_shell, etc).

In this case one can use the built-in Manifold objects, such as SphericalManifold (as you use), PolarManifold, etc.

However, I'm still not sure how to create an internal interface with my own specified function. E.g. something

as simple as creating the domain [-1,1]x[0,1] with the interface $\Gamma$ parameterized by points along the

curve $(0.2 \sin(2\pi y), y)$.

I guess the FunctionManifold class is my best bet, but I'm still trying to understand how this works.

If you have any further suggestions, they would be appreciated.

But thank you again for the helpful comment.

Best regards,

-Sean

[Sean Carney]

Hi all,

Just wanted to follow up on this in case it could help someone else looking to do something similar.

Based on Winni's very helpful response, I was able to successfully create a mesh of the unit square

$[0,1]^2$ that contains the interface described by $\Gamma(x) = (x, 0.5 + 0.1 \sin(2\pi x))$.

Please see the attached--it's not a standalone .cpp file (i.e. not a 'minimal working example'), but, it 

should be able to be copied directly into Step-49 and work right away.

Thanks again--

--Sean