Step-90 TraceFEM tutorial

[vyushut]

Adds a new tutorial step-90 for the trace finite element method (TraceFEM).
Two aspects are of our focus. First, the surface approximation is separated from the discretization of the surface PDE,
e.g. a $Q_2$ discrete level-set and a $Q_1$ solution are possible, but both should be based on the same bulk triangulation.
Second, we make sure that performance of TraceFEM in the parallel implementation corresponds to that of a classical
fitted FEM for a two-dimensional problem.

[tjhei]

Thanks, Vladimir. Can you upload a PDF of the documentation to help with review? You might be able to add it as a file here by dropping it inside the comment field.

[vyushut]

The deal.II Library_ The step-90 tutorial program.pdf

[vyushut]

It seems that the build problems are related to how TrilinosWrappers and linear algebra vectors are used in step-90.cc. Should I use the LA namespace from step-40 and actually follow its code more closely to try to fix the issues?

[tjhei]

Can you fix the CMakeLists.txt to require Trilinos? See step-32. I think the problems will go away then.

[tjhei]

/rebuild

[bangerth]

Here's a bunch of comments about the introduction and results sections.

[bangerth]

Want to state your affiliation and the grant number(s) from NSF?

[vyushut]

@tjhei , could you please provide the appropriate grant number?

[bangerth]

Before you jump right into it in the next sentence, I think this would be a good place to describe what TraceFEM actually is, and what kinds of problems it solves. Perhaps also to reference a paper or two.

In fact, before you describe the method itself, give an outline of the problem you want to solve, and what is difficult about it.

[vyushut]

Here is what I've added:
TraceFEM solves PDEs posed on a possibly evolving $(dim-1)$-dimensional surface employing a fixed uniform background mesh of a $dim$-dimensional domain in which the surface is embedded. Such surface PDEs arise in problems involving material films with complex properties and in other situations in which a non-trivial condition is imposed on either a stationary or a moving interface. Here we consider a steady, complex, non-trivial surface and the prototypical Laplace-Beltrami equation which is a counterpart of the Poisson problem on flat domains.

Being an unfitted method, TraceFEM allows to circumvent the need of remeshing of an evolving surface if it is implicitly given by the zero contour of a level-set function. At the same time, it easily provides with a natural extension of the discrete solution which must be constructed by any time-stepping technique. Certainly, this flexibility comes with the price: one needs to design the nodes and weights for a quadrature customized for each implicit intersection of the zero level-set and the background mesh. Moreover, these intersections may be of arbitrary shape and size manifesting in the so-called ``small cut" problem and requiring addition of a stabilization form that restores well-conditioning of the problem.

[bangerth]

Yes, I think that's great!

In "provides with a natural extension of the discrete solution" I suppose you mean "extension in the direction perpendicular to the surface"? (If so, say "provides a natural", dropping the "with".)

"comes with the price" -> "comes with a price"

[vyushut]

Ahh, I see... I did mean "an extension". I rewrote it to shift the focus as follows:
At the same time, it easily provides with an extension of the discrete solution to a neighborhood of the surface which turns out to be very handy in case of non-stationary interfaces and films.

[bangerth]

Since you compare with "fitted FEM", perhaps that's something that could go into the description of the method I'm asking for above as well: How does TraceFEM compare to other approaches, such as fitted FEM? What do each of these methods do, and how are they different?

[vyushut]

A natural alternative to TraceFEM in solving surface PDEs is the parametric surface finite element method. The latter
method relies on an explicit parametrization of the surface which may be not feasible especially for evolving interfaces with an unknown in advance shape - in this sense, TraceFEM is a technique inspired by the level-set description of interfaces. However, the parametric surface finite element method, when applicable, enjoys many well-known properties of fitted methods on flat domains provided the geometric errors - which are present for both methods - are taken under control.

[bangerth]

What does EOC stand for? And why would you show EOC on the number of DoFs?

[vyushut]

You are right, this is not a "experimental order of convergence". I'll use "Rate" instead. The message of this rate is that it corresponds to the rate of a two-dimensional problem. For a fitted 3D problem we should see something like -3.00. There was another column titled "LevelSet DOFS" with the same 2D rate -2.00 despite the fact that the mesh is three-dimensional. Would you like to see it here too?

[bangerth]

Ah, I see. Perhaps put that interpretation into text below the table to make it clear what you want the interpretation to be.

[bangerth]

Can you explain here what you show in "Stabilization"? The term as shown in the introduction is a bilinear form, not a single number.

[vyushut]

I've changed the title of the column from "Stabilization" to $s_h^{1/2}(u_h)$ and added the following explanation:
...the first order convergence of the stabilization $s_h^{1/2}(u_h)=\sqrt{s_h(u_h, u_h)}$ evaluated at the solution $u_h$.

[bangerth]

Let me know once you have this part too.

[bangerth]

Separately, you will also need to add the necessary lines in the various places in doc/doxygen/tutorial/tutorial.h.in.

[bangerth]

First, thanks a lot for writing this program!

Here is some commentary about the program itself. I made my way down to the start of the main class, but it's late and I'll leave it here for now. I'll come back to it in a few days; go ahead and see whether you can't already address some of these comments.

As a general rule (which we have not always followed in all tutorial programs, but have found correct anyway): We've generally found that it is difficult for people to see the overall flow of a program if comments are interspersed too frequently. As a consequence, we've tried to move comments so that they cover larger blocks of code, rather than just a few lines at a time. I think step-5 is, for instance, a good example of what works well (https://dealii.org/developer/doxygen/deal.II/step_5.html#step_5-CommProg) whereas step-7 is not (https://dealii.org/developer/doxygen/deal.II/step_7.html#step_7-CommProg).

In your program, you often have 1-3 line long comments break apart larger blocks of code. That would be the right strategy if the goal were to just document a piece of code by itself (i.e., a code that is not run through doxygen to generate a website), but it's going to make it hard to read for users and keep the bigger picture in mind. Take a pass through your program and see whether perhaps some of the commentary within a function could be moved to a larger comment up front.

[bangerth]

Use consistent indentation. Throughout the documentation, we try to consistently spell things "DoF" and "DoFs". Perhaps do a search-replace here too.

[vyushut]

Search-replace is done:)

[bangerth]

Let's call this argument IteratorType to make clear it's a type. Do you need this to be a template parameter, though, or would typename DoFHandler<dim>::cell_iterator (or ...::active_cell_iterator) be enough?

[vyushut]

The stabilization should be evaluated on active cells only so typename DoFHandler<dim>::active_cell_iterator is enough.

[tjhei]

First, thanks a lot for writing this program!

Thank you, Wolfgang. Are you okay with the general idea/setup/test problem?

[bangerth]

Yes, sure. I haven't worked my way through the rest of the program (and was hoping that @vyushut could go through the rest and apply some of the general comments first), but I think it's a nice method to solve a nontrivial problem. With a bit more background in the introduction, this should be a nice addition to the tutorial!

[bangerth]

It's ok to leave it as is, but do want to point out that that's a mathematician's approach :-) From an applied perspective, the goal is to solve equations that are relevant because they describe real situations in physics, chemistry, or some other discipline. The equation is what it is, i.e., what came out of a modeling process. Whether or not we like the equation is not a relevant question, and similarly we cannot just choose to add a term to an equation to make the problem easier to solve.

That said, this is a tutorial in which we teach techniques, not solve real applications, and so I'm ok with this statement :-)

[vyushut]

I'm not sure if it helps, but I expanded on why the extra term is included.
We would like to solve the simplest possible problem defined on a surface, namely the Laplace--Beltrami equation, @f{equation*} -\Delta_\Gamma u + c u = f \qquad \text{in }\, \Gamma, @f} where we take $c=1$ for concreteness. We added the term $cu$ to the left-hand side so the problem becomes well-posed in the absence of any boundary; an alternative could be to take $c=0$ but impose the zero mean condition.

[bangerth]

Yes, I think this is great. Thanks for indulging me :-)

[bangerth]

@vyushut Thanks for working through these comments. Leave a comment here (at the top level of the PR) once you've worked through everything and pushed a new version. I will then look through the remainder of the program.

[peterrum]

I like the tutorial 👍 I left some comments.

Maybe @simonsticko could also review the tutorial after you have integrated @bangerth's and my comments.

[peterrum]

Suggested change

	a mesh which is localized near the zero contour line, i.e near the surface, to restore the two-dimensional complexity.
	a mesh along the zero contour line to restore the two-dimensional complexity.

[vyushut]

I wouldn't say "a mesh along the zero contour line" as it indicates that we know in advance which cells need to be created and refined. Instead the mesh is constructed iteratively and a good part of tutorial is devoted to it. I suggest the following:
...a mesh which is localized near the zero contour line of the level set function, i.e near the surface...

[peterrum]

Suggested change

	is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create
	is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create implicitly

[vyushut]

I've rephrased the entire paragraph

[peterrum]

The name is chosen unfortunately!? I guess it is needs_stabilization()!? Also, I am not sure if we need a separate class for the stabilization. You only have a single implementation... Let's keep the cope simple and inline the functions here.

[vyushut]

You are right, there is only one implementation of the stabilization. The method needs_cell_worker() indicates that a bulk cell quadrature must be constructed in addition to the non_matching_fe_values constructed to integrate over implicitly given intersections. Hence the name.
However, there are different TraceFEM stabilizations, e.g. the so called ghost penalty term, which requires fe_values objects ( on faces AND\OR on cells depending on the fe degree) . I think that having a separate class is illuminating and it implements an important logic of separating mesh_workers that require different fe_values from the assembly().

[peterrum]

What is shown? I guess you have created the iso-surface in Paraview?

I think it would make also sense to show 2D results, then you could show the result at the background DoFs. Or alternately, show the results in 3D but with the mesh intersected by plane or so.

[peterrum]

You also need to add the tutorial to https://github.com/dealii/dealii/blob/master/doc/doxygen/tutorial/tutorial.h.in.

Could you also add some references, e.g., to you paper?

[vyushut]

I've added some references.

[simonsticko]

Maybe @simonsticko could also review the tutorial after you have integrated bangerth's and my comments.

Sure. Please tell me if/when you want more comments @vyushut .

[vyushut]

Hello @peterrum . Thank you for the comments! I will them incorporate them and let you know.

[tjhei]

Suggested change

	DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
	constraints.clear();
	constraints.reinit(locally_relevant_dofs);
	DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs(), locally_relevant_dofs);
	constraints.reinit(locally_owned_dofs, locally_relevant_dofs);

[vyushut]

The deal.II Library_ The step-90 tutorial program.pdf