How to extract nodes/vertices at the interface of two domains?

[Behrooz Karami]

Dear all,

I have a mesh composed of two (sub)domains that they share vertices at the interface. Assuming having one Triangulation (and one DoFHandler) for the whole mesh, and in the most simplest case (i.e. no refinement, no hanging nodes, etc) what is a good way to extract such data? To the best of my understanding  there are already a few approaches/methods in Deal.ii used for  extracting similar information (cells and faces).

Thanks very much in advance.

Behrooz Karami

[Jean-Paul Pelteret]

Dear Behrooz,

I guess that the first thing that you need to be able to do is to traverse the mesh until you can get an iterator and index to the cell face that you’re interested in. If the manifold already has some sort of ID associated with it, then you can just do something like

types::manifold_id manifold_of_interest = …;

for (auto cell : dof_handler.active_cell_iterators())

  for (const auto &face : cell->face_iterators())

    if (face->at_boundary()  == false && face->manifold_id() == manifold_of_interest)

    { do something…; }

Otherwise you could maybe get to the manifold face by marking both subdomains that you’ve mentioned with a different material_id and then do something like

types::material_id subdomain_1_mat_id = …;

types::material_id subdomain_2_mat_id = …;

for (auto cell : dof_handler.active_cell_iterators())

  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)

    if (cell->face(f)->at_boundary()  == false && cell->material_id() == subdomain_1_mat_id && cell->neighbor(f)->material_id() == subdomain_1_mat_id)

    { do something…; }

Or you can just use some geometric argument to get to the manifold.

 for (auto cell : dof_handler.active_cell_iterators())

  for (const auto &face : cell->face_iterators())

    if (face->at_boundary()  == false && is_point_on_manifold(face->center()) == true)

    { do something…; }   

Once you have the manifold face that you’re interested in, you can get the vertex positions like this:

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

  const Point<dim> vertex_position = face->vertex(v);

  // OR const Point<dim> vertex_position = cell->face(f)->vertex(v);

If its the solution "at the vertices” that you’re looking for, then you can adapt the code listed in this FAQ entry:

https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-get-the-degree-of-freedom-indices-at-vertices

I hope that this helps you solve your problem!

Best,

Jean-Paul

--
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/5fe1ebd5-ad50-45f3-aa48-aac1d1eddebbn%40googlegroups.com.

[Behrooz Karami]

Dear Jean-Paul,

Thank you very much! That was exactly what I needed to make sure about.

Regards,
Behrooz