Importing nodal BCs and accessing vertices

[Alex Quinlan]

Hi everyone,

I am working on adapting some abaqus jobs for deal.ii.  In the abaqus format, It is typical to specify boundary conditions in terms of (node number, dof number, magnitude).

I have found a working approach for doing this, but I am looking to find a more efficient way as I may end up with millions of nodes and thousands of vertices. Here's an example of how I have been applying nodal displacement BCs:

// define a "node" by the coordinates

    const Point<dim> proot1(0., 0.0209, -0.0101);

  

    // loop through all cells      

    for (const auto &cell : dof_handler.active_cell_iterators())   

    {

         // loop through all of the vertices in a given cell     

         for (unsigned int v = 0; v<cell->n_vertices() ; v++)                {   

              // check if a vertex is near to the node.
              if (cell->vertex(v).distance(proot1) < 0.001)   //              {                                                                      // apply a constraint                                               constraints.add_line(cell->vertex_dof_index(v,0, cell->active_fe_index() ));    

} } }

In an ideal world, I could compute the global dof number based on the node number and local dof number.  Something like:

   global_dof = node_number * 3 + local_dof

   constraints.add_line(global_dof)

However, I understand that this approach isn't implemented in deal.ii because it would be useless after mesh refinement.                                                

Is there a good way to find the global_dof numbers for a vertex/node? 

If not, can I iterate through all of the vertices directly (rather than going through the cell)?  Something like:

for (const auto &vertex : dof_handler.vertex_iterators())

{

         if (vertex.distance(proot1) < 0.001)    

         {

                     constraints.add_line(vertex->vertex_dof_index(0) );                   

}}                               
                                                                                             
 

I'm interested on your thoughts on how to efficiently import nodal boundary conditions from an external file. I couldn't find any other posts on this specific topic, but I apologize if this is a duplicate.                                                         

Thanks,

Alex

[Wolfgang Bangerth]

On 1/13/23 09:55, Alex Quinlan wrote:
>
> Is there a good way to find the global_dof numbers for a vertex/node?

No, and it's not just because of mesh refinement, but also because you may
want to renumber degrees of freedom as implemented in the DoFRenumbering
namespace. For example, if you have a vector displacement problem, you may
want to number DoFs in such a way that the three displacements at each node
are numbered consecutively, or you may want all x-displacements to be numbered
before all y-displacements before all z-displacements. In other words, at
least *in general*, there is no way to tell what the mapping from vertex to
DoF index is.

There is of course also the issue that DoFs live on only vertices only if you
are using linear elements. For quadratic elements, they also live on edge and
face midpoints.

> If not, can I iterate through all of the vertices directly (rather than going
> through the cell)?  Something like:
>
> for (const auto &vertex : dof_handler.vertex_iterators())
> {
>          if (vertex.distance(proot1) < 0.001)
>          {
>                      constraints.add_line(vertex->vertex_dof_index(0) );
> }}

I'm afraid there is no such way, though you can make your life a bit cheaper
if you do something along the lines of

std::vector<bool> vertex_already_handled (triangulation.n_vertices(),
false);

// loop through all cells
for (const auto &cell : dof_handler.active_cell_iterators())
{
// loop through all of the vertices in a given cell
for (unsigned int v = 0; v<cell->n_vertices() ; v++)

if (vertex_already_handled[cell->vertex_index(v)] == false)
{
if (condition)
register constraint;

vertex_already_handled[cell->vertex_index(v)] = true;
}
}

There is of course also the question of whether you have to do the code
snippet above for every imported node, i.e., whether all of the above is
enclosed in a loop of the form

for (unsigned int boundary_node=0; boundary_node<n_nodes; ++...)

It may be cheaper to move that loop into the loop over all cells, but you're
still ending up with a quadratic complexity in the number of cells/boundary
nodes. If that turns out to be too expensive, you need to sort your boundary
nodes into something like an rtree or similar data structure so that looking
up which boundary nodes are close to a cell's vertex becomes substantially
cheaper.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Alex Quinlan]

Thanks Wolfgang.  I appreciate the feedback.  I'll see if I can implement some of the cost-saving tips that you suggested.

[Alex Quinlan]

Hi Wolfgang,

I've come back to this question after discovering some Grid tools.  I am now assessing GridTools::find_closest_vertex and GridTools::find_active_cell_around_point

const Point<3> pt1(0.0, 15.0,  8.20754);
std::vector<bool>  marked_vertices(triangulation.n_vertices(), true);

// Find the vertex closest to the point

unsigned int closest_v1 = GridTools::find_closest_vertex( dof_handler, pt1);

// Mark the vertex for a faster cell search

marked_vertices[closest_v1] = true;

// Find the cell containing the point (or a cell containing, if there are multiple)

auto cell_and_ref_point = GridTools::find_active_cell_around_point(
          dof_handler,    pt1,  marked_vertices, 1.e-4);

//  Search the vertices, then apply constraints

  for (unsigned int v = 0; v<cell_and_ref_point->n_vertices() ; v++)  // might be able to use 'vertex_iterator'
         if (cell_and_ref_point->vertex(v).distance(load1) < 1.e-4)
                constraints.add_line(cell_and_ref_point->vertex_dof_index(v,0, cell_and_ref_point->active_fe_index() ));

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?

Thanks,

Alex

[Wolfgang Bangerth]

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.

[Alex Quinlan]

Thanks, Wolfgang.  I will abide by those guidelines.

[Alex Quinlan]

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

[Luca Heltai]

If you have access to the DG, you can call dh.get_triangulation().n_active_cells(). 

Take a look at the tests subdirectory to find more inspiration on more ways to exploit your information. 

For example, I’d use DoFTools::map_dofs_to_support_points, use that to generate an r tree of indices, and then just ask for the one nearest neighbor (in the tests subdirectory there are a few examples that show how to do this for collections of points). 

Luca

Il giorno 11 apr 2024, alle ore 15:43, Alex Quinlan <alex.r....@gmail.com> ha scritto:



--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1076fa66-e8f5-47a1-bb8c-19cd64e85b2bn%40googlegroups.com.

[Alex Quinlan]

Thanks, Luca.  I will try to implement your suggestions.

Best regards,

Alex