From fa8d3a26cc6b52f3378af5bd7b2f1a14db120037 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 17 Dec 2010 15:54:03 +0000 Subject: [PATCH] Talk about how to extend the program in the direction of slightly crazier domains. git-svn-id: https://svn.dealii.org/trunk@22995 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-38/doc/results.dox | 104 +++++++++++++++++++++++ 1 file changed, 104 insertions(+) diff --git a/deal.II/examples/step-38/doc/results.dox b/deal.II/examples/step-38/doc/results.dox index a3eda5aec0..64b6c89994 100644 --- a/deal.II/examples/step-38/doc/results.dox +++ b/deal.II/examples/step-38/doc/results.dox @@ -35,3 +35,107 @@ Finally, the program produces graphical output that we can visualize. Here is a plot of the results: @image html step-38.solution.png + + + +

Possibilities for extensions

+ +Computing on surfaces only becomes interesting if the surface is more +interesting than just 5/6th of a sphere. To achieve this, deal.II can read +meshes that describe surfaces through the usual GridIn class. Or, in case you +have an analytic description, a simple mesh can sometimes be stretched and +bent into a shape we are interested in. + +Let us consider a relatively simple example: we take the 5/6th sphere we used +before, we stretch it by a factor of 10 in the z-direction, and then we jumble +the x- and y-coordinates a bit. Let's show the computational domain and the +solution first before we go into details of the implementation below: + +@image html step-38.warp-1.png + +@image html step-38.warp-2.png + +The way to produce such a mesh is by using the GridTools::transform +function. It needs a way to transform each individual mesh point to a +different position. Let us here use the following, rather simple function +(remember: stretch in one direction, jumble in the other two): + +@code +template +Point warp (const Point &p) +{ + Point q = p; + q[dim-1] *= 10; + + if (dim >= 2) + q[0] += 2*std::sin(q[dim-1]); + if (dim >= 3) + q[1] += 2*std::cos(q[dim-1]); + + return q; +} +@endcode + +If we followed the LaplaceBeltrami::make_mesh function, we would +extract the 5/6th spherical surface mesh as before, warp it into the shape we +want, and refine as often as necessary. This is not quite as simple as we'd +like here, though: refining requires that we have an appropriate manifold +object attached to the triangulation that describes where new vertices of the +mesh should be located upon refinement. I'm sure it's possible to describe +this manifold in a not-too-complicated way by simply undoing the +transformation above (yielding the spherical surface again), finding the +location of a new point on the sphere, and then re-warping the result. But I'm +a lazy person, and since doing this is not really the point here, let's just +make our lives a bit easier: we'll extract the 5/6th sphere, refine it as +often as necessary, get rid of the object that describes the manifold since we +now no longer need it, and then finally warp the mesh. With the function +above, this would look as follows: + +@code +template +void LaplaceBeltrami::make_mesh () +{ + HyperBallBoundary boundary_description; + Triangulation volume_mesh; + GridGenerator::half_hyper_ball(volume_mesh); + + volume_mesh.set_boundary (1, boundary_description); + volume_mesh.set_boundary (0, boundary_description); + volume_mesh.refine_global (6); + + static HyperBallBoundary surface_description; + triangulation.set_boundary (1, surface_description); + triangulation.set_boundary (0, surface_description); + + std::set boundary_ids; + boundary_ids.insert(0); + + GridTools::extract_boundary_mesh (volume_mesh, triangulation, + boundary_ids); + triangulation.set_boundary (1); /* ** */ + triangulation.set_boundary (0); /* ** */ + GridTools::transform (&warp, triangulation); /* ** */ + + std::cout << "Surface mesh has " << triangulation.n_active_cells() + << " cells." + << std::endl; +} +@endcode + +Note that the only addition has been the three lines marked with asterisks. It +is worth pointing out one other thing here, though: because we un-attach the +manifold description from the surface mesh, whenever we use a mapping object +in the rest of the program, it has no curves boundary description to go on any +more. Rather, it will have to use the implicit, StraightBoundary class that is +used on all parts of the boundary not explicitly assigned a different +mannifold object. Consequently, whether we use MappingQ(2), MappingQ(15) or +MappingQ1, each cell of our mesh will be mapped using a bilinear +approximation. + +All these drawbacks aside, the resulting pictures are still pretty. The only +other differences to what's in step-38 is that we changed the right hand side +to $f(\mathbf x)=\sin x_3$ and the boundary values (through the +Solution class) to $u(\mathbf x)|_{\partial\Omega}=\cos x_3$. Of +course, we now non longer know the exact solution, so the computation of the +error at the end of LaplaceBeltrami::run will yield a meaningless +number. -- 2.39.5