From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Wed, 10 Sep 2014 00:46:14 +0000 (-0500)
Subject: Add a section discussing consequences.
X-Git-Tag: v8.2.0-rc1~155^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=77382897996843300ab2053675106b092f42bff7;p=dealii.git

Add a section discussing consequences.

This pertains in particular to the question of getting discontinuous
geometries when using manifolds and a discussion of what that means,
if it is important, and if it could be fixed.
---

diff --git a/examples/step-53/doc/results.dox b/examples/step-53/doc/results.dox
index 6f7dc62860..beaaf746e8 100644
--- a/examples/step-53/doc/results.dox
+++ b/examples/step-53/doc/results.dox
@@ -120,3 +120,163 @@ values, we just go down this list and set these elevations in the
 their original zero value). This vector is then output using DataOut as
 usual and can be visualized as shown above.
 
+
+<h3>Issues with adaptively refined meshes generated this way</h3>
+
+If you zoomed in on the mesh shown above and looked closely enough, you would
+find that at hanging nodes, the two small edges connecting to the hanging
+nodes are not in exactly the same location as the large edge of the
+neighboring cell. This can be shown more clearly by using a different surface
+description in which we enlarge the vertical topography to enhance the effect
+(courtesy of Alexander Grayver):
+
+<img src="http://www.dealii.org/images/steps/developer/step-53.smooth-geometry.png" alt="">
+
+So what is happening here? Partly, this is only a result of visualization, but
+there is an underlying real cause as well:
+
+<ul>
+  <li>When you visualize a mesh using any of the common visualization
+  programs, what they really show you is just a set of edges that are plotted
+  as straight lines in three-dimensional space. This is so because almost all
+  data file formats for visualizing data only describe hexahedral cells as a
+  collection of eight vertices in 3d space, and do not allow to any more
+  complicated descriptions. (This is the main reason why
+  DataOut::build_patches() takes an argument that can be set to something
+  larger than one.) These linear edges may be the edges of the cell you do
+  actual computations on, or they may not, depending on what kind of mapping
+  you use when you do your integrations using FEValues. By default, of course,
+  FEValues uses a linear mapping (i.e., an object of class MappingQ1) and in
+  that case a 3d cell is indeed described exclusively by its 8 vertices and
+  the volume it fills is a trilinear interpolation between these points,
+  resulting in linear edges. But, you could also have used tri-quadratic,
+  tri-cubic, or even higher order mappings and in these cases the volume of
+  each cell will be bounded by quadratic, cubic or higher order polynomial
+  curves. Yet, you only get to see these with linear edges in the
+  visualization program because, as mentioned, file formats do not allow to
+  describe the real geometry of cells.
+
+  <li>That said, let us for simplicity assume that you are indeed using a
+  trilinear mapping, then the image shown above is a faithful representation
+  of the cells on which you form your integrals. In this case, indeed the
+  small cells at a hanging nodes do not, in general, snugly fit against the
+  large cell but leave a gap or may intersect the larger cell. Why is this?
+  Because when the triangulation needs a new vertex on an edge it wants to
+  refine, it asks the manifold description where this new vertex is supposed
+  to be, and the manifold description duly returns such a point by (in the
+  case of a geometry derived from ChartManifold) pulling the adjacent points
+  of the line back to the reference domain, averaging their locations, and
+  pushing forward this new location to the real domain. But this new location
+  is not usually along a straight line (in real space) between the adjacent
+  vertices and consequently the two small straight lines forming the refined
+  edge do not lie exactly on the one large straight line forming the unrefined
+  side of the hanging node.
+</ul>
+
+The situation is slightly more complicated if you use a higher order mapping
+using the MappingQ class, but not fundamentally different. Let's take a
+quadratic mapping for the moment (nothing fundamental changes with even higher
+order mappings). Then you need to imagine each edge of the cells you integrate
+on as a quadratic curve despite the fact that you will never actually see it
+plotted that way by a visualization program. But imagine it that way for a
+second. So which quadratic curve does MappingQ take? It is the quadratic curve
+that goes through the two vertices at the end of the edge as well as a point
+in the middle that it queries from the manifold. In the case of the long edge
+on the unrefined side, that's of course exactly the location of the hanging
+node, so the quadratic curve describing the long edge does go through the
+hanging node, unlike in the case of the linear mapping. But the two small
+edges are also quadratic curves; for example, the left small edge will go
+through the left vertex of the long edge and the hanging node, plus a point it
+queries halfway in between from the manifold. Because, as before, the point
+the manifold returns halfway along the left small edge is rarely exactly on
+the quadratic curve describing the long edge, the quadratic short edge will
+typically not coincide with the left half of the quadratic long edge, and the
+same is true for the right short edge. In other words, again, the geometries
+of the large cell and its smaller neighbors at hanging nodes do not touch
+snuggly.
+
+This all begs two questions: first, does it matter, and second, could this be
+fixed. Let us discuss these in the following:
+
+<ul>
+  <li>Does it matter? It is almost certainly true that this depends on the
+  equation you are solving. For example, it is known that solving the Euler
+  equations of gas dynamics on complex geometries requires highly accurate
+  boundary descriptions to ensure convergence of quantities that are measure
+  the flow close to the boundary. On the other hand, equations with elliptic
+  components (e.g., the Laplace or Stokes equations) are typically rather
+  forgiving of these issues: one does quadrature anyway to approximate
+  integrals, and further approximating the geometry may not do as much harm as
+  one could fear given that the volume of the overlaps or gaps at every
+  hanging node is only ${\cal O}(h^d)$ even with a linear mapping and ${\cal
+  O}(h^{d+p-1})$ for a mapping of degree $p$. (You can see this by considering
+  that in 2d the gap/overlap is a triangle with base $h$ and height ${\cal
+  O}(h)$; in 3d, it is a pyramid-like structure with base area $h^2$ and
+  height ${\cal O}(h)$. Similar considerations apply for higher order mappings
+  where the height of the gaps/overlaps is ${\cal O}(h^p)$.) In other words,
+  if you use a linear mapping with linear elements, the error in the volume
+  you integrate over is already at the same level as the integration error
+  using the usual Gauss quadrature. Of course, for higher order elements one
+  would have to choose matching mapping objects.
+
+  Another point of view on why it is probably not worth worrying too much
+  about the issue is that there is certainly no narrative in the community of
+  numerical analysts that these issues are a major concern one needs to watch
+  out for when using complex geometries. If it does not seem to be discussed
+  often among practitioners, if ever at all, then it is at least not something
+  people have identified as a common problem.
+
+  This issue is not dissimilar to having hanging nodes at curved boundaries
+  where the geometry description of the boundary typically pulls a hanging
+  node onto the boundary whereas the large edge remains straight, making the
+  adjacent small and large cells not match each other. Although this behavior
+  existed in deal.II since its beginning, 15 years before manifold
+  descriptions became available, it did not ever come up in mailing list
+  discussions or conversations with colleagues.
+
+  <li>Could it be fixed? In principle, yes, but it's a complicated
+  issue. Let's assume for the moment that we would only ever use the MappingQ1
+  class, i.e., linear mappings. In that case, whenever the triangulation class
+  requires a new vertex along an edge that would become a hanging node, it
+  would just take the mean value of the adjacent vertices <i>in real
+  space</i>, i.e., without asking the manifold description. This way, the
+  point lies on the long straight edge and the two short straight edges would
+  match the one long edge. Only when all adjacent cells have been refined and
+  the point is no longer a hanging node would we replace its coordinates by
+  coordinates we get by a manifold. This may be awkward to implement, but it
+  would certainly be possible.
+
+  The more complicated issue arises because people may want to use a higher
+  order MappingQ object. In that case, the Triangulation class may freely
+  choose the location of the hanging node (because the quadratic curve for the
+  long edge can be chosen in such a way that it goes through the hanging node)
+  but the MappingQ class, when determining the location of mid-edge points
+  must make sure that if the edge is one half of a long edge of a neighboring
+  coarser cell, then the midpoint cannot be obtained from the manifold but
+  must be chosen along the long quadratic edge. For cubic (and all other odd)
+  mappings, the matter is again a bit complicated because one typically
+  arranges the cubic edge to go through points 1/3 and 2/3 along the edge, and
+  thus necessarily through the hanging node, but this could probably be worked
+  out. In any case, even then, there are two problems with this:
+
+  - When refining the triangulation, the Triangulation class can not know what
+    mapping will be used. In fact it is not uncommon for a triangulation to be
+    used differently in different contexts within the same program. If the
+    mapping used determines whether we can freely choose a point or not, how,
+    then, should the triangulation locate new vertices?
+  - Mappings are purely local constructs: they only work on a cell in
+    isolation, and this is one of the important features of the finite element
+    method. Having to ask whether one of the vertices of an edge is a hanging
+    node requires querying the neighborhood of a cell; furthermore, such a
+    query does not just involve the 6 face neighbors of a cell in 3d, but may
+    require traversing a possibly very large number of other cells that
+    connect to an edge. Even if it can be done, one still needs to do
+    different things depending on how the neighborhood looks like, producing
+    code that is likely very complex, hard to maintain, and possibly slow.
+
+  Consequently, at least for the moment, none of these ideas are
+  implemented. This leads to the undesirable consequence of discontinuous
+  geometries, but, as discussed above, the effects of this do not appear to
+  pose problem in actual practice.
+
+</ul>