... 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);
Magdalena schrieb am Freitag, 16. Dezember 2022 um 11:09:40 UTC+1: > 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 > > <http://Marching%20Cube%20Algorithm>). > > 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 > > marco....@gmail.com schrieb am Donnerstag, 15. Dezember 2022 um 13:26:21 > UTC+1: > >> 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 >> >> Il giorno giovedì 15 dicembre 2022 alle 05:55:51 UTC+1 Praveen C ha >> scritto: >> >>> 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 >> >> -- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/531b3eee-faeb-4410-80be-ab4da3547914n%40googlegroups.com.
