Transfinite Interpolation

[Narendra Nanal]

Hello All,

I am trying to run benchmark case of 2D flow around a cylinder. Here are the details- 

http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark1_re20.html.

My question is how to use Transfinite Interpolation class to generate a mesh with a smooth transition. Thank you. 

Regards,

Narendra 

[David Wells]

Hi Narendra

We should probably put something like this in step-49. Here is a bit of code that shows how to use the TFI manifold class:

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>

int main()
{
  using namespace dealii;

  Triangulation<2> tria;
  GridGenerator::hyper_cube_with_cylindrical_hole (tria, 0.05, 0.1);
  const types::manifold_id tfi_id = 1;
  for (const auto cell : tria.active_cell_iterators())
    {
      // set all non-boundary manifold ids to the new TFI manifold id.
      cell->set_manifold_id(tfi_id);
      for (unsigned int face_n = 0; face_n < GeometryInfo<2>::faces_per_cell;
           ++face_n)
        {
          if (!cell->face(face_n)->at_boundary())
            cell->face(face_n)->set_manifold_id(tfi_id);
        }
    }
  TransfiniteInterpolationManifold<2> inner_manifold;
  inner_manifold.initialize(tria);
  tria.set_manifold(1, inner_manifold);

  tria.refine_global(2);

  GridOut go;
  std::ofstream out("grid.eps");
  go.write_eps(tria, out);
}

This creates a nice blend between the inner circle and the outer square. From here you can do something like what is done in step-49 to generate the rest of the grid from Cartesian pieces.

If you are using 8.5 instead of 9.0 you will need to make all Manifold classes (including putting your own PolarManifold on the circle) before the Triangulation.

I hope this helps!

Thanks,
David Wells

[Martin Kronbichler]

Hi Narenda and David,

Just to add on top of what David Wells said, Luca Heltai and I once had the idea of showing this manifold in a new tutorial program (where we combine it with MappingFEField to show the performance impact). Now step-60 does contain some MappingFEField, but on a different topic, so I think we should simply create such a program soon. I will open an issue for that.

Best,
Martin

--
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+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[David Wells]

Hi Martin and Narenda,

I was curious to see how hard it would be to put together a nice grid for the DFG benchmark problem. Its rather involved: we should look into a way to extend the library so that coming up with high quality grids like this one is too hard. Anyway, here is what I came up with for the offset cylinder:

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>

#include <iostream>
#include <fstream>

int main()
{
  using namespace dealii;

  // set up the bulk triangulation
  Triangulation<2> bulk_triangulation;
  GridGenerator::subdivided_hyper_rectangle(bulk_triangulation,
                                            {22u, 4u},
                                            Point<2>(0.0, 0.0),
                                            Point<2>(2.2, 0.41));
  std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
  Tensor<1, 2> cylinder_triangulation_offset;
  for (const auto cell : bulk_triangulation.active_cell_iterators())
    {
      if ((cell->center() - Point<2>(0.2, 0.2)).norm() < 0.15)
        cells_to_remove.insert(cell);

      if (cylinder_triangulation_offset == Point<2>())
        {
          for (unsigned int vertex_n = 0; vertex_n < GeometryInfo<2>::vertices_per_cell; ++vertex_n)
            if (cell->vertex(vertex_n) == Point<2>())
              {
                // skip two cells in the bottom left corner
                cylinder_triangulation_offset = 2.0*(cell->vertex(3) - Point<2>());
                break;
              }
        }
    }
  Triangulation<2> result_1;
  GridGenerator::create_triangulation_with_removed_cells(bulk_triangulation,
                                                         cells_to_remove,
                                                         result_1);

  // set up the cylinder triangulation
  Triangulation<2> cylinder_triangulation;
  GridGenerator::hyper_cube_with_cylindrical_hole (cylinder_triangulation, 0.05, 0.41/4.0);
  GridTools::shift(cylinder_triangulation_offset, cylinder_triangulation);
  // dumb hack
  for (const auto cell : cylinder_triangulation.active_cell_iterators())
    cell->set_material_id(2);

  // merge them together
  auto minimal_line_length = [](const Triangulation<2> &tria) -> double
    {
      double min_line_length = 1000.0;

      for (const auto cell : tria.active_cell_iterators())
        {

          min_line_length = std::min(min_line_length,
                                     (cell->vertex(0) - cell->vertex(1)).norm());
          min_line_length = std::min(min_line_length,
                                     (cell->vertex(0) - cell->vertex(2)).norm());
          min_line_length = std::min(min_line_length,
                                     (cell->vertex(1) - cell->vertex(3)).norm());
          min_line_length = std::min(min_line_length,
                                     (cell->vertex(2) - cell->vertex(3)).norm());
        }
      return min_line_length;
    };

  // the cylindrical triangulation might not match the Cartesian grid: as a
  // result the vertices might not be lined up. Get around this by deleting
  // duplicated vertices with a very low numerical tolerance.
  const double tolerance = std::min(minimal_line_length(result_1),
                                    minimal_line_length(cylinder_triangulation))/2.0;


  Triangulation<2> result_2;
  GridGenerator::merge_triangulations(result_1,
                                      cylinder_triangulation,
                                      result_2,
                                      tolerance);

  const types::manifold_id tfi_id = 1;

  const types::manifold_id polar_id = 0;
  for (const auto cell : result_2.active_cell_iterators())

    {
      // set all non-boundary manifold ids to the new TFI manifold id.

      if (cell->material_id() == 2)
        {

          cell->set_manifold_id(tfi_id);
          for (unsigned int face_n = 0; face_n < GeometryInfo<2>::faces_per_cell;
               ++face_n)
            {

              if (cell->face(face_n)->at_boundary())
                cell->face(face_n)->set_manifold_id(polar_id);
              else
                cell->face(face_n)->set_manifold_id(tfi_id);
            }
        }
    }

  PolarManifold<2> polar_manifold(Point<2>(0.2, 0.2));
  result_2.set_manifold(polar_id, polar_manifold);
  TransfiniteInterpolationManifold<2> inner_manifold;
  inner_manifold.initialize(result_2);
  result_2.set_manifold(tfi_id, inner_manifold);

  std::vector<Point<2> *> inner_pointers;
  for (const auto cell : result_2.active_cell_iterators())
    {

      for (unsigned int face_n = 0; face_n < GeometryInfo<2>::faces_per_cell;
           ++face_n)
        {

          if (cell->face(face_n)->manifold_id() == polar_id)
            {
              inner_pointers.push_back(&cell->face(face_n)->vertex(0));
              inner_pointers.push_back(&cell->face(face_n)->vertex(1));
            }
        }
    }
  // de-duplicate
  std::sort(inner_pointers.begin(), inner_pointers.end());
  inner_pointers.erase(std::unique(inner_pointers.begin(),
                                   inner_pointers.end()),
                       inner_pointers.end());

  // find the current center...
  Point<2> center;
  for (const Point<2> *const ptr : inner_pointers)
    center += *ptr/double(inner_pointers.size());
  std::cout << "computed center: "
            << center
            << '\n';

  // and recenter at (0.2, 0.2)
  for (Point<2> *const ptr : inner_pointers)
    *ptr += Point<2>(0.2, 0.2) - center;

  Point<2> center2;
  for (const Point<2> * const ptr : inner_pointers)
    center2 += *ptr/double(inner_pointers.size());
  std::cout << "recomputed center: "
            << center2
            << '\n';

  GridOut go;
  std::ofstream out("grid.eps");

  result_2.refine_global(1);
  go.write_eps(result_2, out);
}

I attached a picture of the output. This script relies on a patch that I need to add to deal.II that permits merging triangulations with a tolerance (so it won't work yet on 9.0, but hopefully it will soon!).

Thanks,

David Wells

[Wolfgang Bangerth]

On 05/26/2018 07:56 AM, David Wells wrote:
>
> I was curious to see how hard it would be to put together a nice grid for the
> DFG benchmark problem.

After the other patch is merged, why don't you put that into the library as
well? It's a common benchmark, and we have GridGenerator::cheese(). We may as
well as GridGenerator::DFG_flow_around_cylinder_geometry() or similar.

Cheers
W.

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