Sharp-interface test

[peterrum]

This PR adds a test to compute the force vector based on the sharp-interface method:

With a bit of cosmetic in deal.II, the code will look like this:

// gather-evaluate curvature
cell->get_dof_values(curvature_solution, buffer);
phi_curvature.evaluate(cell, unit_points, buffer, EvaluationFlags::values);

for (int i = 0; i < dim; ++i)
  {
    // gather-evaluate normal i
    cell->get_dof_values(normal_vector_field.block(i), buffer);
    phi_normal_force.evaluate(cell, unit_points, buffer, EvaluationFlags::values);

    // loop over all quadrature points and perform surface integral
    for (unsigned int q = 0; q < n_q_points; ++q)
      phi_normal_force.set_value(phi_normal_force.value(q) * phi_curvature.value(q) * JxW[q]);

    // integrate-scatter force i
    phi_normal_force.integrate(cell, unit_points, buffer, EvaluationFlags::values);
    cell->distribute_local_to_global(buffer, force_vector.block(i));
  }

which matches exactly the formula :100: I guess we could skip the inner loop (but not sure about this due to BlockVector).

I still need to clean up the code and sort the quadrature points according to cells so that we can exploit the full potential of FEPointValues. Anyway, cool feature!

Currently, the surface mesh is not moving. It is generated with the pre-knowledge of initial condition of the level-set field (which is a lot comparing that I did not have a clue how to proceed yesterday evening).

Test 1: computation of curvature and normal based on LS

... with black being quantities on the background mesh and red quantities on the surface mesh

Test 2: computation of curvature and normal based on geometric properties of surface mesh

Test 3: identify volume enclosed by surface mesh

Test 4: rising bubble

Solve rising-bubble test case with approach 2.

Test 5: move surface mesh by velocity of background mesh

[peterrum]

I have added a second test. It focuses how to extract normal and curvature at the surface quadrature points only by relying on geometrical information and without using the level-set field. For this purpose, I am using MappingFEField ((absolute) Eulerian description of positions of the support points of elements), I am computing the normal via FEValues::normal_vector() and the curvature FEValues::get_function_gradients(). ~I actually think that I can use for the operations FEEvaluation.~ Nope: dim==spacedim!

Since MappingFEField relies on a vector of positions of the support points, it is probably quite simple to update the geometry description of the surface mesh (by "moving" the points).

The following algorithms could be derived: evaluate velocity at the surface quadrature points on the background mesh, move the surface mesh (by updating MappingFEField), and the quadrature point, normal vector n, curvature kappa and quadrature weight JxW (or the product of these n kappa JxW) to the background mesh and test them there.

[peterrum]

@mschreter The failing test is related to updates in deal.II. There is a major clean up of the ReferenceCell/ReferenceCells namespace. I hope it is over soon and the interface stabilizes. But, not sure when that will happen. Let's not update deal.II for now!

[peterrum]

@mschreter Thanks for the comments here and offline! I hope we will be able to use parts of the implementations in MeltPoolDG in the future!

[peterrum]

@mschreter @kronbichler Thanks for the review! You both have approved the additions. I am still a bit hesitant in merging in master, since I still cannot guarantee that the implementation 100% correct. But if you think, we can nevertheless merge (the next PR wouldn't contain 4k loc additions).

[peterrum]

@mschreter @kronbichler I have squashed all commits. If you feel it is ready for master, feel free to merge.

[peterrum]

Let me merge this PR. I would like to work on a small project that was in the pipline for quite some (and is not SIM related) for which I can use the new "level-set solver".