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
On Thursday, May 24, 2018 at 8:43:38 AM UTC-4, Martin Kronbichler wrote:
>
> 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
>
> On 23.05.2018 23:32, David Wells wrote:
>
> 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
>
>
>
> On Monday, May 21, 2018 at 7:06:19 PM UTC-4, Narendra Nanal wrote:
>>
>> 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
>>
> --
> 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 [email protected] <javascript:>.
> For more options, visit https://groups.google.com/d/optout.
>
>
>
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
