Hi Wolfgang,
I've gotten an R-tree subroutine set up for applying my boundary conditions
based on rtree_8h.html
<https://www.dealii.org/current/doxygen/deal.II/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
<https://www.dealii.org/current/doxygen/deal.II/tria__description_8cc.html#ab695339a227bd369fcea767631247cff>
.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.
