[cig-commits] [commit] master: Use a manifold to generate better meshes for the shell. (f050b9b)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/aspect

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/23371c522d510d3edddde05d16731fd7039d510b...eef6299b61e1952be22819a80f6a6e4ee344b677

>---------------------------------------------------------------

commit f050b9b216eebe0d3a3544ee286c73412f363e59
Author: Wolfgang Bangerth <bangerth at math.tamu.edu>
Date:   Fri Oct 31 13:37:29 2014 -0500

    Use a manifold to generate better meshes for the shell.
    
    This only works with a sufficiently new deal.II, of course.


>---------------------------------------------------------------

f050b9b216eebe0d3a3544ee286c73412f363e59
 doc/modules/changes.h                    | 13 ++++++++++++-
 source/geometry_model/spherical_shell.cc | 32 +++++++++++++++++++++++++++++++-
 2 files changed, 43 insertions(+), 2 deletions(-)

diff --git a/doc/modules/changes.h b/doc/modules/changes.h
index 2258038..d1c691f 100644
--- a/doc/modules/changes.h
+++ b/doc/modules/changes.h
@@ -6,6 +6,18 @@
  *
  *
  * <ol>
+ * <li> New: ASPECT now uses meshes for the spherical shell where each node is
+ * placed on concentric shells. Previously, upon mesh refinement, new vertices
+ * were placed at an averaged location of the vertices of the mother cell,
+ * leading to a certain degree of distortion of cells. It also led to a
+ * situation where not all cells are equally shaped due to this
+ * distortion. The new mesh, in contrast, is completely
+ * symmetric. Addtionally, a higher order mapping to represent curved faces is
+ * now also used for the interior cells using the manifold description of the
+ * shell/sphere.
+ * <br>
+ * (Wolfgang Bangerth, Timo Heister, 2015/01/16)
+ *
  * <li> New: GeometryModel::has_curved_elements() will allow for optimizations
  * if all boundaries (and manifold description) are straight.
  * <br>
@@ -39,7 +51,6 @@
  * output time step will produce output. Old checkpoint files will continue
  * to work, only with a possible short gap in visualization output right
  * after restart.
- * 
  * <br>
  * (Rene Gassmoeller, 2014/12/03)
  *
diff --git a/source/geometry_model/spherical_shell.cc b/source/geometry_model/spherical_shell.cc
index 5195f64..0c1f8e8 100644
--- a/source/geometry_model/spherical_shell.cc
+++ b/source/geometry_model/spherical_shell.cc
@@ -23,7 +23,10 @@
 
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria_boundary_lib.h>
-
+#if (DEAL_II_VERSION_MAJOR>=9) || \
+  ((DEAL_II_VERSION_MAJOR==8) && (DEAL_II_VERSION_MINOR >= 2))
+#  include <deal.II/grid/manifold_lib.h>
+#endif
 
 namespace aspect
 {
@@ -76,6 +79,33 @@ namespace aspect
           Assert (false, ExcInternalError());
         }
 
+#if ((DEAL_II_VERSION_MAJOR>=9) ||				\
+     ((DEAL_II_VERSION_MAJOR==8) && (DEAL_II_VERSION_MINOR >= 2)))
+
+      // if deal.II is sufficiently new, then use a manifold
+      // description for all cells. use manifold_id 2 in order not to
+      // step on the boundary indicators used below
+      static const SphericalManifold<dim> spherical_manifold;
+      coarse_grid.set_manifold (2, spherical_manifold);
+
+      for (typename Triangulation<dim>::active_cell_iterator
+	     cell = coarse_grid.begin_active();
+	   cell != coarse_grid.end(); ++cell)
+	cell->set_all_manifold_ids (2);
+
+      // clear the manifold id from objects for which we have boundary
+      // objects (and need boundary objects because at the time of
+      // writing, only boundary objects provide normal vectors)
+      for (typename Triangulation<dim>::active_cell_iterator
+	     cell = coarse_grid.begin_active();
+	   cell != coarse_grid.end(); ++cell)
+	for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+	  if (cell->at_boundary(f))
+	    cell->face(f)->set_all_manifold_ids (numbers::invalid_manifold_id);
+#endif
+      
+      // attach a boundary object to the inner
+      // and outer boundaries
       static const HyperShellBoundary<dim> boundary_shell;
       coarse_grid.set_boundary (0, boundary_shell);
       coarse_grid.set_boundary (1, boundary_shell);