Added GridTools::get_subdomain_association(Tria, CellId).

[marcfehling]

Blocked by #11364 and #11365.

The function returns the subdomain of a cell represented by a CellId on a parallel::distributed::Triangulation. We ask the p4est oracle for the corresponding owner.

For all other triangulation types, we simply request the subdomain id of a cell iterator object created on the local object from a CellId. @peterrum -- This works for fullydistributed too, right?

I still need to come up with a meaningful test for the new functionality, which is why I open this as a draft. But I would already welcome your thoughts and opinions!

[peterrum]

For what do you such a function?

@peterrum -- This works for fullydistributed too, right?

no

[marcfehling]

We would like to create a map of cells which live on different triangulations that may be partitioned differently. This function will be used to determine sending and receiving processes for data transfer between those cells. To be more specific, between cells on a volume and a surface mesh.

[peterrum]

We would like to create a map of cells which live on different triangulations that may be partitioned differently.

In this case, I add a Do not merge label. I don't think this is a good approach - going through p4est cell by cell seem to be too expensive. We should discuss alternative approaches. I am doing something similar here in https://github.com/dealii/dealii/blob/master/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h.

[bangerth]

The p4est case is actually quite cheap. You're creating the p4est_quadrant object, which is O(# of levels), and then the look-up via p4est_comm_owner is O(log_2 (# of processors)). In particular, it doesn't iterate over all cells to find which the owner is. In other words, the overall complexity of this function is O(log_2 (# of processors) + log(N)) because # of levels = O(log N).

The point of the function (which should probably be mentioned in its documentation) is that with the function one can also ask for the owner of cells identified by CellId that don't actually exist on the current processor. So for the f::d::T, this is going to be expensive. I'd implement this with an Assert(false, ExcNotImplemented()) for the moment.

[bangerth]

This looks good to me, but it needs to be re-based once #11364 is merged.

@peterrum We need this kind of functionality and it's cheap to query for p::d::T without communication. At some point we probably want this to also work with p::f::T. Are you ok with merging this (once rebased)? Do you have an idea how at some point one would implement this also for p::f::T?

[peterrum]

I had a look at the implementation of the p4est function (https://github.com/cburstedde/p4est/blob/ea25d0221dce486dc40ea224a394b06c9f016a6b/src/p4est_communication.c#L675-L742). Apparently, they save the global start point of each process on each process. Via binary search they can find the rank of a cell.

This only works because they exploit the structure of a space-filling curve. This is not something we can do in the case other partitioning approaches like done in p:f:T. In this case, we would need to do a global communication per cell if the function is called for each cell instead for a batch.

Furthermore, I am not happy with Triangulation getting a new function (we should rather try to decrease the number of functions of Triangulation). For me this is more a GridTools function.

As a consequence, I would rather suggest to have a function with the following interface:

namespace GridTools
{
  template<int dim, int spacedim>
  std::vector<unsigned int>
  get_owning_subdomain_for_cells(const Triangulation<dim, spacedim>&  tria, 
                                 const std::vector<CellId> &cell_id);
}

where one could batch the query. For p:d:T, the function would use the implementation as suggested in this PR. For p:f:T, there would be no implementation for now. For p:s:T, the artificial cells would be removed as a preparation. Furthermore, the documentation could say that it is highly recommended not to call this function for a single cell.

Are you ok with merging this (once rebased)?

As stated above, I am not happy with this PR in the current state. But I won't block.

[bangerth]

@marcfehling originally had the function externally, but it clearly requires internal knowledge of the triangulation in question to work -- so he moved it into the triangulation class, which I believe is the right place for it to be.

I'd be ok with changing the interface so that it take a whole vector of CellId objects. That would certainly work for the use case we have in mind. Would that alleviate your concerns?

[peterrum]

I'd be ok with changing the interface so that it take a whole vector of CellId objects. That would certainly work for the use case we have in mind. Would that alleviate your concerns?

That would be better.

But, I am still not convinced that the serial Triangulation is the right place for such a function. For me, the function belongs more to the (re)partitioning algorithm which is (currently) part of the derived classes.

[kronbichler]

But, I am still not convinced that the serial Triangulation is the right place for such a function. For me, the function belongs more to the (re)partitioning algorithm which is (currently) part of the derived classes.

I would second this point - this is a very specific use case that seems to break with our previous tradition of keeping features genuine to a parallel/distributed triangulation out of the base class. To give an example, the much more widely used function get_communicator is not part of dealii::Triangulation, and this causes some significant effort in our code base. A quick grep for dynamic_cast<const parallel::TriangulationBase as well as the related dealii::parallel::TriangulationBase gives 70 hits, of which maybe 60 are simply to get a communicator for various downstream consumers. So if we add this function to the base class, I do not see why we have refrained from allowing access to the communicator from the base class without dynamic casts.

[tjhei]

I do not see why we have refrained from allowing access to the communicator from the base class without dynamic casts.

I am happy to accept a PR that adds this.

[bangerth]

Let's not hold this PR up because we've made design mistakes in other contexts :-)

I'm not quite sure what the alternative to the function here would be, though. Would you want to introduce it as a virtual function parallel::Triangulation? I get that that would make sense if the implementation in the base class can not produce a reasonable result -- like in the case of get_mpi_communicator(). But here, the base class can generate a reasonable answer, and it would avoid having to try whether a particular triangulation is derived from parallel::Triangulation or not when deciding whether we can call the function.

Maybe it's worthwhile explaining where we want to use this function: We'd like to write something of the sort

  template <typename VolumeTriangulationType, typename SurfaceTriangulationType>
  XXX
  extract_surface_mesh (const VolumeTriangulationType &volume_mesh, 
                                        SurfaceTriangulationType &surface_mesh);

The returned object of type XXX will need to contain a map from volume faces to surface cell and back, in some shape of form. This is not difficult if both triangulation types are sequential: Then we could just store regular iterators. But in the parallel case, we can't identify cells with iterators if they are not stored on the current processor -- we have to use CellId for this task, and then we will have to set up a data structure that allows us to communicate, for example, data extracted from the faces of the volume mesh to the owners of the corresponding surface cell or the other way around. These kinds of maps and communication strategies require finding out who owns what.

The point is that it would be quite nice if we could set all of this up in some sort of generic way. Being able to query the owner of a cell identified by CellId is a key operation we have to figure out somehow, and it is an operation that has a natural implementation for a sequential triangulation: The answer is "process 0", which is what the base class should implement.

[kronbichler]

Thanks for sharing the application background, this makes the discussion easier, even though I have to admit I could probably have extracted that from #11366 (comment) already.

What I think @peterrum is concerned about is the generality of this approach, and alternatives on a higher level, not on where to place this function. Given

dealii/include/deal.II/grid/tria_accessor.h

Lines 3485 to 3486 in 97bbe3d

	types::global_cell_index
	global_active_cell_index() const;

and with a ConsensusAlgorithm, one can identify this information for any kind of triangulation. One merely needs an intermediate step with bounding boxes to do it in a fully unstructured way for cases where p4est cannot help in terms of a complicated coarse grid. As a matter of fact, we have done something very similar to the coupling between a volume and surface mesh for a fluid-structure-interaction application with non-matching grids along the fluid-structure interface already, coded mostly by @peterrum regarding the communication aspect, as a result of many discussions with me and @nfehn. So I would really like to join the efforts here and aim to insert the most general approach in the library. Or a small set of different variants, in case it turns out the flexible setup were not as suited in terms of HPC for example (I am not suggesting the two are in conflict at this stage, just a random example). For this, we would need to define what exactly the surface mesh and volume mesh set as requirements and how we would like the communication pattern to look like. I could not find that specification in #10037 so I suggest we take that effort now to save a lot of work down the road.

[marcfehling]

I do not have a general opinion on where to add this function.

@bangerth's point to make it a class member is that we need access to internal data structures. While this is true, we can also get access to those via getter functions that are already implemented (most probably). His other argument to make it a member of the dealii::Triangulation class is that there is a valid answer to this query for any type of Triangulation. This allows for a similar interface among all triangulation classes, but this functionality will most probably never be used in the sequential context.

There are enough arguments for and against it.

As far as I can tell, we allow the GridTools::partition_triangulation_* to take any type of Triangulation (excluding p::d). So personally, I would also add this particular functionality for any type of Triangulation, including the sequential case. I do not know whether to place this function in the GridTools namespace or to make it a class member though.

---

I do however have an opinion on the generality of this approach.

For every different Triangulation type in the parallel namespace, I could think of a different implementation that is optimized for their use case: The p::s case just needs to return the corresponding entry of the true_subdomain_ids_of_cells member. In the p::d case we can look up the owner with the binary search in comm_find_owner which does not involve any communication between other processors. The p::f case needs the Consensus approach. I would prefer if we had an optimized algorithm for every triangulation type.

However, the interface for all functions should look the same. @peterrum suggested to supply a container of CellId objects, which I would be fine with.

Since the Consensus approach should work for any kind of parallel Triangulation, I would suggest that we should implement it for either the parallel::TriangulationBase or parallel::DistributedTriangulationBase class, and add specializations for the p::s and p::d cases, e.g. by overwriting them. The p::f class and others to come would then inherit the base implementation.

What do you think?

[bangerth]

@kronbichler We continue to not make progress and are holding up @marcfehling :-(

I understand that if you have non-matching grids, that one needs a more general algorithm that will involve bounding boxes. But we're not interested in this, and in fact have much more information available than just a general cell's location: The CellId object stores the coarse mesh cell and its child-grandchild-... path through the forest. One can write substantially more efficient algorithms for this case, and that's exactly what we're trying to do here -- as evidenced by the fact that with the p4est implementation, this lookup can be done with out communication and with log P complexity. For the p::f::T case, my understanding is that partitioning will be based exclusively on the coarse cell index, so we know that the optimal implementation will be equivalent to one global-reduction with communication with essentially trivial local look-up on every process. That, too, is more efficient than what one can conceive with non-matching meshes.

What is your objection to having a function

  std::vector<unsigned int>
  get_owning_subdomain_for_cells(const std::vector<CellId> &cell_id);

in the interface of Triangulation? I know that you will need a different function for what you have in mind, but these two cases are clearly different, have different complexities, and are therefore consequently best addressed with two different functions.

[nfehn]

I have gone through all the comments of this PR, but I have to admit that it is really difficult to extract the problem description. @bangerth, you argued that you prefer simpler algorithms than some generic routines: What I did not understand so far is where you get the relevant information from to construct the maps if the routines are not "fully generic"? How does the part of your algorithm look like that cares about geometry, location, neighborhood?

@bangerth You have also motivated this PR from the point of view of complexities. For me, complexity is only part of the problem, because identifying and removing bottlenecks (talking about absolute runtimes) is more important. It would be helpful to get insights why you have concerns that a generic algorithm will not satisfy your needs. In my opinion one should first argue or give evidence that the generic routines suggested by @kronbichler form a bottleneck relative to other grid/triangulation functions that one has to do anyway, because two different functions (which are again specialized for all the dealii::triangulations if I understood the discussion correctly) is something that one would like to avoid from a code design, maintenance, documentation, usability perspective. I wonder a bit why the opinions are so divided regarding this question. Isn't the discussion also stuck because the topic is so difficult that one always has to think about the different Triangulation types, that p4est is part of this discussion, etc.? In my opinion, this should be a clear hint to rethink the design.

[kronbichler]

Sorry for the delay.

The main problem I see with this algorithm is that it only realistically works with two types of triangulations, the serial one (where you do not need it) and for the parallel::distributed::Triangulation with p4est. But it is hard or impossible to adapt to fully distributed triangulations or any kind of triangulation with a large coarse mesh, where one needs a different approach without the function proposed here. This is why I would like to not expose it in central classes and rather as a free function in GridTools or in a derived class, which does not fool users into expecting something we cannot realistically deliver. I realize that we would not be able to use p4est information in GridTools without somehow hacking p::d::Tria, but in my opinion the alternative implementation would not be that much more expensive in the grand scheme, so I would prefer it.

As always, I do not feel that strongly that I would straight out oppose to this PR and the last thing I would want to do is block your progress, but I wanted to raise a flag that a different interface might be preferable in terms of later extensions. That said, we would still need to discuss the different requirements on coupling algorithms and as things look like we won't find a one-way-fits-all solution at least in the short run, so a CellId based solution might provide some value. But for that broader discussion #10037 is the right place.

[marcfehling]

I think one of the problems related to this discussion is that we do not have the entire solution for coupling bulk and surface processes in front of us, and still work towards the final picture.

I would thus suggest to place this function in the GridTools namespace for now. It will be one independent function that can easily be modified. The alternative would be deep-rooting it in the Triangulation class hierarchy as a member function, which would make this harder to maintain for the moment.

If we decide that we do not need this function, we could simply scrap it later on. If we have a plan how to acquire this information on many or all types of Triangulation classes, we can still turn it into a member function and deprecate the free function.

I will adjust this PR accordingly.

[peterrum]

/rebuild

[peterrum]

Looks good. Just a few minor points.

[peterrum]

Suggested change

	// the most general type of triangulation is the serial one. here, all
	// subdomain information is directly available
	for (const auto &cell_id : cell_ids)
	{
	subdomain_ids.push_back(
	triangulation.create_cell_iterator(cell_id)->subdomain_id());
	}
	return std::vector<types::subdomain_id>(cell_ids.size(), 0);

[marcfehling]

I would leave this as is.

If, at some point, we decide to handle the subdomain_id for serial triangulations differently (e.g., choose invalid_subdomain_id instead of 0), we need to remember to also change it here.

[marcfehling]

If there are no further objections, I will squash all commits and release this PR from its draft state.

[peterrum]

Sounds good! Thanks!

[marcfehling]

I've squashed all commits. Two comments on the current state:

• I've renamed the function to get_subdomain_association as there is already a similar one in the same namespace so we just overloaded it
• there are only explicit instantiations of this function for dim >= 2 as there is no p::d::Triangulation for dim == 1

[marcfehling]

ping

[luca-heltai]

As far as I can see, all comments have been addressed. I'll merge. :)