Hi Wolfgang,
I've gotten an R-tree subroutine set up for applying my boundary conditions based on
rtree_8h.html. It seems to work, but I have a couple difficulties that I've highlighted in red. Hopefully, others might find use in my adaptation of the example and in my questions.
My subroutine takes the following arguments: the DOF handler, the constraints object, and a vector of a custom nBC structure for holding the nodal boundary conditions that I've read from the file.
template <int dim, int spacedim>
void nodal_bcs_rtree( std::vector<nBC> & nBCs,
dealii::DoFHandler<dim,spacedim> & dof_handler,
dealii::AffineConstraints<double> & constraints )
{
namespace bgi = boost::geometry::index;
// -------- Create the R-tree -----------
std::vector<BoundingBox<3>> all_boxes(100000); // TODO Need a way to set the correct size for this
std::vector<dealii::TriaActiveIterator<dealii::DoFCellAccessor<3, 3, false> >> cell_list(100000);
unsigned int i = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
all_boxes[i] = cell->bounding_box(); // prepare the bounding boxes to be packed into the r-tree
cell_list[i++] = cell; // create a corresponding vector of cells accessors
}
const auto idx_tree = pack_rtree_of_indices(all_boxes);
// -------- Search the R-tree ----------------
const double nrad = 1e-5; // Tolerance on vertex position TODO make relative to cell size
for (nBC bc : nBCs)
{
std::vector<unsigned int> closest_cell_idx;
idx_tree.query(bgi::nearest(bc.coords,1), std::back_inserter(closest_cell_idx)); // bc.coords is type Point<3>
if (closest_cell_idx[0] < i-1) // Check that the cell index number will correspond to an actual cell in cell_list
{
const auto cell = cell_list[closest_cell_idx[0]];
for (unsigned int v = 0; v<cell->n_vertices() ; v++) {
if (cell->vertex(v).distance(bc.coords) < nrad) {
~~~Apply constraints ~~~
}
}
}
}
}
----------------
1) I'm having trouble getting the number of cells, and then sizing my vectors appropriately. From what I can see, the DoFHandler won't provide the number of cells, only the number of DoFs. In the short term, I've just thrown in the number 100,000 as the vector size. In the example ( https://www.dealii.org/current/doxygen/deal.II/rtree_8h.html ), the line:
std::vector<BoundingBox<2>> all_boxes(tria.n_locally_owned_active_cells());
gets used to set the vector size. From within my subroutine, I don't have access to the triangulation object. Am I missing a way to get "n_locally_owned_active_cells" from the DoFHandler? Or should I just pass this value in as an argument to my subroutine?
2) Does my approach of creating vectors to hold the bounding boxes and the cell pointers make sense? Is there a better way to get access to the cell once I've determined the BoundingBox index from the R-tree?
3) Somewhat related to question 1, is there a way to check that the 'closest_cell' was actually found? I had an issue with MPI where the target point 'bc.coords' was outside the bounding box of some processes, resulting in some errors. I set this (closest_cell_idx[0] < i-1) to check that an actual BoundingBox was found.
Thanks for your suggestion on the R-tree. I'd welcome any critiques of my implementation.
Cheers,
Alex
On Wednesday, November 15, 2023 at 8:35:49 AM UTC-5 Alex Quinlan wrote:
Thanks, Wolfgang. I will abide by those guidelines.
On Tuesday, November 14, 2023 at 10:15:29 PM UTC-5 Wolfgang Bangerth wrote:
On 11/14/23 16:25, Alex Quinlan wrote:
>
> I'm curious what your thoughts are on this approach. I imagine it could have
> an advantage in some situations, depending on the size of the mesh and the
> number of constraints to be added.
>
> I have not done any speed testing on this yet, tho. Do you think it would be
> looking into? Or do you see some fatal flaw with this approach?
Alex:
the general approach most professional programmers will ascribe to is to write
a version of the code that is intelligible and easy to maintain. Only then do
you worry about speed. If the code in question is fast enough (say at most of
few percent up to 20% of the program's run time -- as measured with a class
such as TimerOutput), then it's not worth worrying about it.
So the questions you're asking are premature. Make the code work what it is
you want it to do, and then you can think about its performance.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.colostate.edu/~bangerth/
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1076fa66-e8f5-47a1-bb8c-19cd64e85b2bn%40googlegroups.com.
