Finding cut cells

[Praveen C]

Dear all

I am looking for available methods to find cells cut by some given curve.

I found such a case in step-60 and made this example

https://bitbucket.org/cpraveen/deal_ii/src/master/embedded_curve/

As shown in the figures, some cells are not identified though they are cut by the curve. This is a limitation which is already explained in step-60.

Are there other methods available in deal.II for this purpose ? Or do you know of any external library than can be used for this, and which can be called from a deal.II code ?

I am assuming that the curve is given as piecewise polynomials. I only need 2-D case.

Thanks
praveen

[Marco Feder]

Hi, 

I think Step 85 does this by using a (global) level set function. 

Otherwise, you can use RTrees of bounding boxes with https://www.dealii.org/developer/doxygen/deal.II/classGridTools_1_1Cache.html#aca2782d6e93b5a0033c046b57904c67f combined with boost::geometry::index to identify intersections, in case where you don't want to use a level-set and you just work with the mesh itself. The following snippet does this:

namespace bgi = boost::geometry::index;

const auto &tree =

space_cache -> get_locally_owned_cell_bounding_boxes_rtree();

// Bounding boxes of the embedded grid

const auto &embedded_tree = embedded_cache -> get_cell_bounding_boxes_rtree();

for (const auto &[embedded_box, embedded_cell] : embedded_tree)

{

for (const auto &[space_box, space_cell] :

tree | bgi::adaptors::queried(bgi::intersects(embedded_box)))

{

// do what you need here

}

}

Best,

Marco

[Magdalena]

Dear Praveen, 

if you have an implicit description of your curve, let's say a level-set like DoF-vector,

VectorType level_set_vector;

 where a value "contour_value" represents your surface, a helpful tool could be the Marching Cube Algorithm (https://dealii.org/developer/doxygen/deal.II/classGridTools_1_1MarchingCubeAlgorithm.html#ab29b812b67e5aa185c0f5966163f40fa). 

The latter is set up via 

GridTools::MarchingCubeAlgorithm<dim, VectorType> mc(mapping,

                                                     dof_handler.get_fe());

Then, you could loop over the cells in your triangulation, and check whether the processed cell is cut or not: 

std::vector<const typename Triangulation<dim, dim>::cell_iterator> cut_cells

for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_locally_owned())
{

// determine points and cells of aux surface triangulation
std::vector<Point<dim>> surface_vertices;
std::vector<CellData<dim == 1 ? 1 : dim - 1>> surface_cells;

if (dim > 1)
     mc.process_cell(cell, level_set_vector, contour_value, surface_vertices, surface_cells);

// no surface cells are written

else
    mc.process_cell(cell, level_set_vector, contour_value, surface_vertices);

// append cell if it is cut by the interface 

if (surface_vertices.size() == 0)
   cut_cells.emplace_back(cell);

}

}

I hope this matches for your purpose.

Best,
Magdalena

[Magdalena]

... Sorry, there is a mistake in one of the last code lines. It should be:

// append cell if it is cut by the interface 

if (surface_vertices.size() > 0)
   cut_cells.emplace_back(cell);