Finding vertices on boundaries

[Shahab Golshan]

Dear all,

I want to find all the vertices located on boundaries. It seems the function at_boundary() does not work for vertices (points). Do you have any suggestions how to obtain all the vertices located on the boundaries?
Best,

Shahab

[Bruno Blais]

Dear Shahab,

I think the best solution (I might be wrong) is to loop through the boundary faces and acquire the DOF indices that lie on a boundary face. You could then store them in an std::map or a similar structure that prevents data duplication.

You can create a vector and store the dof indices of a face in a following fashion that I believe should work.

Assuming you already have a vector such as:

std::vector<types::global_dof_index> face_dofs;





In the loops within the cell then the faces, assuming you are at a face face and that you have a FiniteElement fe, you can do:

face_dofs.resize(fe.dofs_per_face);
face->get_dof_indices(face_dofs, cell->active_fe_index());

This will gather all the dofs indices of that given face. If that face is a boundary face it means you have gathered all of its dofs location.

I hope that points you in the right direction and helps you!

Best

Bruno

[David Wells]

Hi Shahab,

You are right that vertices themselves don't store any information - this is an implementation detail of deal.II (points do not store any topological information in a Triangulation, unlike lines, quadrilaterals, or hexahedra).

Just to check: when you write 'vertices', are you looking for information about the degrees of freedom on the boundary, or are you interested in the getting a list of the Point<spacedim> objects that are on the boundary? What Bruno describes will do the former. If you want to do the later then something like this will work:

Triangulation<dim, spacedim> triangulation;

// do something to set up the triangulation

std::set<Point<spacedim>> boundary_vertices;

for (const auto &face : triangulation.active_face_iterators())

  if (face->at_boundary())

    for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)

      boundary_vertices.insert(face->vertex(v));

You could also store Point<spacedim> *> objects if you want to manipulate the points once you have all of them.

Does this answer your question? Please post again if you still need help figuring this out!

   

Best,

David

--
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/389defa8-9b6b-467e-a2f4-3d77abee31f0%40googlegroups.com.

[Bruno Blais]

I think David solution is more along the lines of Shahab original question (at least from my understanding) :).

To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[Shahab Golshan]

Hi David,

Yes, it was very helpful and solved the problem. Thank you very much.

Best,

Shahab

[Shahab Golshan]

I have also another question on the function find_cells_adjacent_to_vertex.
As input argument, it needs the vertex index, not the vertex as a point. How can I find the index of a vertex (to use in this function) from the vertex location (point).

Best,

Shahab

[Wolfgang Bangerth]

On 5/1/20 5:07 PM, Shahab Golshan wrote:
> I have also another question on the function */find_cells_adjacent_to_vertex/*.

> As input argument, it needs the vertex index, not the vertex as a point. How
> can I find the index of a vertex (to use in this function) from the vertex
> location (point).

cell->vertex_index(v)

is your friend, where v=0...2^d :-)

Best
W.


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

[Doug Shi-Dong]

Hello everyone,

Sorry to revive this thread, but I have a similar question that follows up his question. In the case of a distributed grid, there might be vertices/DoFs on a corner (or even edges in 3D) of a cell. 

For example, see attached figure (at the trailing edge of the airfoil in the square, in the yellow square). Domain 0 is not aware that one of its corner is on a boundary. To aggravate the issue, the DoFHandler with FE_Q distributes the DoFs such that Domain 0 owns that corner DoF.

It seems like I might have to loop over the neighbors of every cell to make sure that account for all the corners/edges DoFs that might be on a boundary, but that seems inefficient.

Anyone knows of a better way within deal.II to deal with this?

Best regards,

Doug

[Wolfgang Bangerth]

On 10/21/20 1:05 PM, Doug Shi-Dong wrote:
>
> It seems like I might have to loop over the neighbors of every cell to make
> sure that account for all the corners/edges DoFs that might be on a boundary,
> but that seems inefficient.
>
> Anyone knows of a better way within deal.II to deal with this?

Each vertex of the locally owned part of the domain is surrounded by either
locally owned cells or ghost cells. As a consequence, if you loop over all
locally owned cells and all ghost cells, you will find that vertex as one
vertex on a boundary face.

You are right that it is difficult to do this by starting at one vertex. The
easier approach is to first loop over all locally-owned + ghost cells, loop
over their faces, see whether a face is on the boundary, and then collect all
vertices that you are interested in in a std::set or similar. Then, later on
when you are curious about one specific vertex, you only need to look whether
it's in the std::set.

[Doug Shi-Dong]

Thanks, yes, that should work well for 2D. I actually just thought of a corner case (pun intended) in 3D. It is possible for a cell's vertex to be on a boundary, without any of its neighbour actually have a boundary face.

Imagine a 3D cube with a hollow inner cube. The cells directly diagonal to the corner cells will have a vertex on the boundary, but none of its neigbor actually lie on the boundary.

Any advice on how to approach this since the current processor would not find it within its ghost cells?

Doug

[Wolfgang Bangerth]

On 10/21/20 7:46 PM, Doug Shi-Dong wrote:
> Thanks, yes, that should work well for 2D. I actually just thought of a corner
> case (pun intended) in 3D. It is possible for a cell's vertex to be on a
> boundary, without any of its neighbour actually have a boundary face.

Yes, precisely: None of the face neighbors have to be on the boundary. But at
least one of the vertex neighbors have to be -- which is why the ghost layer
consists of all of the vertex neighbors of locally owned cells.

You can create the same situation in 2d, by the way: Start with a single
vertex and put 5 or more quadrilateral cells around it. Then put another layer
of cells around those. Then remove one of the cells that are adjacent to the
original vertex and look at the cell opposite from it.

> Any advice on how to approach this since the current processor would not find
> it within its ghost cells?

Yes, it will -- see above.

[Doug]

Thanks for the clarification! I initially thought ghost neighbours would only include face neighbours. That solves the problem.