[cig-commits] [commit] master: Fixed initial conditions for box with origin (ff44215)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/c1465eb80ac07139659936dd80414c93d876e518...1a6b7735396a24c6b2d09b93f5acbcbe86bc2e64

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

commit ff44215f5e39cbb7faf27df1391767e88d59d303
Author: Ryan Grove <rgrove at clemson.edu>
Date:   Tue Dec 9 14:10:38 2014 -0500

    Fixed initial conditions for box with origin


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

ff44215f5e39cbb7faf27df1391767e88d59d303
 include/aspect/geometry_model/box.h |  8 +++++++-
 source/geometry_model/box.cc        |  6 ++++++
 source/initial_conditions/box.cc    | 14 +++++++-------
 3 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/include/aspect/geometry_model/box.h b/include/aspect/geometry_model/box.h
index 22b52c5..e8027ab 100644
--- a/include/aspect/geometry_model/box.h
+++ b/include/aspect/geometry_model/box.h
@@ -46,12 +46,18 @@ namespace aspect
         void create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const;
 
         /**
-         * Return a point that denotes the upper right corner of the box
+         * Return a point that denotes the size of the box in each dimension of the
          * domain.
          */
         Point<dim> get_extents () const;
 
         /**
+         * Return a point that denotes the lower left corner of the box
+         * domain.
+         */
+        Point<dim> get_origin () const;
+
+        /**
          * Return the typical length scale one would expect of features in
          * this geometry, assuming realistic parameters.
          *
diff --git a/source/geometry_model/box.cc b/source/geometry_model/box.cc
index 9b0815d..31daf91 100644
--- a/source/geometry_model/box.cc
+++ b/source/geometry_model/box.cc
@@ -141,6 +141,12 @@ namespace aspect
       return extents;
     }
 
+    template <int dim>
+    Point<dim>
+    Box<dim>::get_origin () const
+    {
+      return box_origin;
+    }
 
     template <int dim>
     double
diff --git a/source/initial_conditions/box.cc b/source/initial_conditions/box.cc
index c3206fd..61d3382 100644
--- a/source/initial_conditions/box.cc
+++ b/source/initial_conditions/box.cc
@@ -42,7 +42,7 @@ namespace aspect
 
       double perturbation = 1;
       for (unsigned int d=0; d<dim; ++d)
-        perturbation *= std::sin(numbers::PI*position[d]/geometry->get_extents()[d]);
+        perturbation *= std::sin(numbers::PI*(position[d]-geometry->get_origin()[d])/geometry->get_extents()[d]);
       return 1 + perturbation/10;
     }
 
@@ -62,8 +62,8 @@ namespace aspect
       Point<dim> temporary1, temporary2;
       for (int d=0; d<dim; ++d)
         {
-          temporary1[d]=geometry->get_extents()[d]*0.625;
-          temporary2[d]=geometry->get_extents()[d]*0.375;
+          temporary1[d]=geometry->get_extents()[d]*0.625+geometry->get_origin()[d];
+          temporary2[d]=geometry->get_extents()[d]*0.375+geometry->get_origin()[d];
         }
 
       return 1+(1/exp(position.distance(temporary2)) - 1/exp(position.distance(temporary1)));
@@ -89,8 +89,8 @@ namespace aspect
       if (center[1] < ratio)
         ratio = center[1];
 
-      double zx = (position[0] - center[0])/ratio;
-      double zy = (position[1] - center[1])/ratio;
+      double zx = (position[0] - geometry->get_origin()[0] - center[0])/ratio;
+      double zy = (position[1] - geometry->get_origin()[1] - center[1])/ratio;
       double x = zx;
       double y = zy;
 
@@ -161,13 +161,13 @@ namespace aspect
             }
           else if (inclusion_gradient == "linear")
             {
-              perturbation = ((radius - position.distance (center)) / radius) * (inclusion_temperature - ambient_temperature);
+              perturbation = ((radius - position.distance(center)) / radius) * (inclusion_temperature - ambient_temperature);
             }
           else if (inclusion_gradient == "constant")
             {
               perturbation = inclusion_temperature - ambient_temperature;
             }
-          if (position.distance (center) > radius)
+          if (position.distance(center) > radius)
             perturbation = 0;
         }