Update point_is_in_domain function

[anne-glerum]

We have use cases where we apply traction boundary conditions and mesh deformation. However, the initial lithostatic pressure traction plugin uses the function point_is_in_domain of the geometry model to check whether the reference point where the lithostatic profile is calculated lies in the domain. At the moment, the function cannot deal with initial topography and throws an error message. It also neglects initial mesh deformation. However, for chunk and box geometries, we know what initial topography is applied. And we can also check whether initial mesh deformation is applied. This PR:

• includes max initial topography at the queried point in the domain extents for box and chunk geometries
• checks whether mesh deformation is applied in the box model. If so, it will throw when initial mesh deformation is applied, or if time-dependent mesh deformation is applied after timestep 0

The problem is that when the lithostatic profile is computed, initialization of the simulator has not finished, and therefore we cannot ask for the Mesh Deformation Handler. My original problem is therefore not solved. I do consider this PR an improvement by itself and would appreciate input on how to also fix my problem. I could read in the required mesh deformation information from the input file, but that seems less elegant.

The initial lithostatic pressure plugin actually does not need to know whether the representative point lies within the domain boundaries in all directions, only in the horizontal directions, as the vertical coordinate is ignored. The point_is_in_domain function could therefore become a point_is_within_horizontal_extents_of_domain function.
But, the point_is_in_domain function is by now used in one other place than for the initial lithostatic pressure plugin, namely in initial_temperature/random_gaussian_perturbation.cc. Therefore, I could also do the checking of the horizontal extents within the initial lithostatic pressure plugin, as this plugin requests a lot of information from the geometry models already.

For all pull requests:

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[gassmoeller]

Hi Anne, I will review this later (have to finish AGU poster first :-)), but here is an idea for your original problem:

• The reason the geometry model cannot determine a point after mesh deformation is that we currently rely on global information about the mesh to determine this (extent, origin, original shape).
• However there is another way, and that is to check all cells of the triangulation whether they contain the point to test. This is expensive if you do it for many points (it scales with N_pointstotest * N_cells), but this plugin only needs to test a single point so it shouldnt be too slow.
• So instead of throwing an assertion if mesh deformation is enabled, we could fall back to the slower method to determine if a point is inside the domain. There is even a deal.II function for it (find_active_cell_around_point, https://dealii.org/developer/doxygen/deal.II/namespaceGridTools.html#a2e10aeb1c8e76110a84b6945eac3aaf0). I think the deal.II function throws an assertion if the point is not in the domain so we may have to do a try ... catch to see if it finds the point. Also this approach has the disadvantage that before the mesh is deformed, we of course dont know how it will be deformed in the future.

Is it worth trying to implement the fall back method in case mesh deformation is enabled? Or is your question specifically about the case where initial mesh deformation will deform the mesh a lot?

[anne-glerum]

Ah good idea! I'll have a look if using find_active_cell_around_point can cover all the possible scenarios, thanks!

[gassmoeller]

I took a closer look at the code, but I have to say it is quite hard to understand its intent. I have two question:

• for the box geometry: is the only change in this PR that now point_is_in_domain can be used in timestep 0 (after simulator initialization) with mesh deformation unless the ascii data mesh deformation is active?
• for the chunk geometry: Why can you remove the assert there, and not add the same conditions as for the box?

If the solution to expand the functionality of the functions to also check deformed meshes is working, I think that would be what I prefer to the current approach in this PR.

[anne-glerum]

There are two things that I wanted to change:

1. instead of throwing when initial topography is applied, use the initial topography at the given point to augment the domain extents. This way we can check whether the point is in the domain for simulations that have initial topography.
2. Expand the cases for which we can determine whether the point is in the domain for simulations with mesh deformation.

Concerning point 1:
For the box geometry, we can ask what the maximum topography is at a certain horizontal position and just add it to the vertical extent. Then we can check as before whether the point falls within the domain extents.
The chunk geometry already takes care of topography when 'pulling back' a point, so the Assert that checks whether initial topography is applied, could be removed.

Concerning point 2:
Time-dependent mesh deformation is not applied during t0, so this is a case where the point_is_in_domain function is still useful. Initial mesh deformation without further time-dependent mesh deformation could be treated the same way as initial topography. Therefore I tried to find out what mesh deformation is applied. For my specific use-case, this didn't help though, as the initial_lithostatic_pressure plugin calls this function before the simulator is initialized.

However, now I took your suggestion and when there is mesh deformation of any kind and at any time, I use the GridTools function. This means all situations are covered (I still need to test it though).
For the chunk, I can replace the remaining Assert that checks for mesh deformation with the same find_active_cell_around_point bit.

[anne-glerum]

I've now also added the find_active_cell_around_point call to chunk.cc. Two new tests check whether the point_is_in_domain function works with initial topography for the box and the chunk. And I've adapted an existing test to test the function over time with mesh deformation. All tests fail on my laptop though, so I have to see what the tester says.

[bangerth]

I am not sure if I like the call to the GridTools function.

[anne-glerum]

Ah yes, of course. I've done something similar to point_values.cc now, is that what you're looking for?

[anne-glerum]

The new tests run now (but fail as I haven't copied over the test results yet). dealii-master seems broken, 18 unrelated gmg tests fail.

[gassmoeller]

Yes, I think this is the right approach, thanks for the changes. Just some smallish comments. Also please add a changelog, its worth noting that this works now (at least for some geometries) with mesh deformation.

[anne-glerum]

Hi @gassmoeller, I've added a changelog, unified the comments and moved part of the functionality into a Utilities function. Anything else? The failing tester is only for gmg tests on dealii main.

[gassmoeller]

Thank you. Looks good to go now. 👍