[cig-commits] [commit] master: Use higher order mapping in interior, correct manifold for shell (546de43)

[cig_noreply at geodynamics.org]

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

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

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

commit 546de433c8368efd051df52703de9c069a4a90ec
Author: Timo Heister <timo.heister at gmail.com>
Date:   Sun Jan 11 20:32:19 2015 -0500

    Use higher order mapping in interior, correct manifold for shell
    
    - set straight boundary for quarter/half shell walls
    - only set boundary objects as needed for spherical shell
    - remove deal.II version check
    - use higher order mapping also in the interior


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

546de433c8368efd051df52703de9c069a4a90ec
 source/geometry_model/spherical_shell.cc | 46 ++++++++++++++++----------------
 source/simulator/core.cc                 |  8 ++++--
 2 files changed, 29 insertions(+), 25 deletions(-)

diff --git a/source/geometry_model/spherical_shell.cc b/source/geometry_model/spherical_shell.cc
index 0c1f8e8..b374d93 100644
--- a/source/geometry_model/spherical_shell.cc
+++ b/source/geometry_model/spherical_shell.cc
@@ -23,10 +23,7 @@
 
 #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
+#include <deal.II/grid/manifold_lib.h>
 
 namespace aspect
 {
@@ -79,33 +76,36 @@ 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
+      // Use a manifold description for all cells. use manifold_id 99 in order
+      // not to step on the boundary indicators used below
       static const SphericalManifold<dim> spherical_manifold;
-      coarse_grid.set_manifold (2, spherical_manifold);
+      coarse_grid.set_manifold (99, 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);
+           cell = coarse_grid.begin_active();
+           cell != coarse_grid.end(); ++cell)
+        cell->set_all_manifold_ids (99);
 
       // 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
+           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);
+
+      // deal.II wants boundary objects even for the straight boundaries
+      // when using manifolds in the interior:
+      static const StraightBoundary<dim> straight_boundary;
+      std::set<types::boundary_id> ids = get_used_boundary_indicators();
+      for (std::set<types::boundary_id>::iterator it = ids.begin();
+           it!=ids.end(); ++it)
+        if (*it > 1)
+          coarse_grid.set_boundary (*it, straight_boundary);
+
+      // attach boundary objects to the curved boundaries:
       static const HyperShellBoundary<dim> boundary_shell;
       coarse_grid.set_boundary (0, boundary_shell);
       coarse_grid.set_boundary (1, boundary_shell);
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 48daf95..bcbe660 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -155,12 +155,16 @@ namespace aspect
 
     //Fourth order mapping doesn't really make sense for free surface
     //calculations (we disable curved boundaries) or we we only have straight
-    //boundaries
+    //boundaries. So we either pick MappingQ(4,true) or MappingQ(1,false)
     mapping (
       (parameters.free_surface_enabled
        ||
        geometry_model->has_curved_elements() == false
-      )?1:4),
+      )?1:4,
+      (parameters.free_surface_enabled
+       ||
+       geometry_model->has_curved_elements() == false
+      )?false:true),
   
     // define the finite element. obviously, what we do here needs
     // to match the data we provide in the Introspection class