[cig-commits] [commit] master: move more inclusion stuff out (c4f333d)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/c883fc069e0dcc93cce45e1d52c84c05c9c4385b...1f4ab820b0191d9a4058b78c91b629bc94c437ec

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

commit c4f333da37f6940f01cb3bfe98baa38f8cb8c702
Author: Timo Heister <timo.heister at gmail.com>
Date:   Tue May 20 11:38:26 2014 -0400

    move more inclusion stuff out


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

c4f333da37f6940f01cb3bfe98baa38f8cb8c702
 benchmark/inclusion/inclusion.cc   | 87 +++++++++++++++++++++++++++++++++++++-
 source/postprocess/duretz_et_al.cc | 37 +---------------
 2 files changed, 87 insertions(+), 37 deletions(-)

diff --git a/benchmark/inclusion/inclusion.cc b/benchmark/inclusion/inclusion.cc
index 33435d9..9a849a9 100644
--- a/benchmark/inclusion/inclusion.cc
+++ b/benchmark/inclusion/inclusion.cc
@@ -8,6 +8,79 @@ namespace aspect
   {
     using namespace dealii;
 
+
+    namespace AnalyticSolutions
+    {
+      // based on http://geodynamics.org/hg/cs/AMR/Discontinuous_Stokes with permission
+      void _Inclusion(double pos[2], double r_inclusion, double eta, double *vx, double *vy, double *p)
+      {
+        const double min_eta = 1.0;
+        const double max_eta = eta;
+        const double epsilon = 1; //strain rate
+        const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
+        std::complex<double> phi, psi, dphi;
+        const double offset[2]= {1.0, 1.0};
+        double r2_inclusion = r_inclusion * r_inclusion;
+
+        double x = pos[0]-offset[0];
+        double y = pos[1]-offset[1];
+        double r2 = x*x+y*y;
+
+        std::complex<double> z(x,y);
+        if (r2<r2_inclusion)
+          {
+            //inside the inclusion
+            phi=0;
+            dphi=0;
+            psi=-4*epsilon*(max_eta*min_eta/(min_eta+max_eta))*z;
+          }
+        else
+          {
+            //outside the inclusion
+            phi=-2*epsilon*A*r2_inclusion/z;
+            dphi=-phi/z;
+            psi=-2*epsilon*(min_eta*z+A*r2_inclusion*r2_inclusion/(z*z*z));
+          }
+        double visc = (r2<r2_inclusion)? max_eta : 1.0;
+        std::complex<double> v = (phi - z*conj(dphi) - conj(psi))/(2.0*visc);
+        *vx=v.real();
+        *vy=v.imag();
+        *p=-2*epsilon*dphi.real();
+      }
+    }
+
+    namespace VelocityBoundaryConditions
+    {
+      namespace DuretzEtAl
+      {
+        template <int dim>
+        Inclusion<dim>::Inclusion ()
+          :
+          eta_B (1e3)
+        {}
+
+
+
+        template <int dim>
+        Tensor<1,dim>
+        Inclusion<dim>::
+        boundary_velocity (const Point<dim> &p) const
+        {
+          Assert (dim == 2, ExcNotImplemented());
+
+          double pos[2]= {p(0),p(1)};
+
+          Tensor<1,dim> velocity;
+          double pressure;
+          aspect::DuretzEtAl::AnalyticSolutions::_Inclusion
+          (pos,0.2,eta_B, &velocity[0], &velocity[1], &pressure);
+
+          return velocity;
+        }
+      }
+    }
+
+
       /**
        * A material model that describes the "Pure shear/Inclusion" benchmark
        * of the paper cited in the documentation of the DuretzEtAl namespace.
@@ -382,10 +455,22 @@ namespace aspect
 {
   namespace MaterialModel
   {
- ASPECT_REGISTER_MATERIAL_MODEL(Inclusion,
+     ASPECT_REGISTER_MATERIAL_MODEL(Inclusion,
                                      "Inclusion",
                                      "A material model that corresponds to the 'Inclusion' benchmark "
                                      "defined in Duretz et al., G-Cubed, 2011.")
  
   }
+  namespace VelocityBoundaryConditions
+  {
+    namespace DuretzEtAl
+    {
+      ASPECT_REGISTER_VELOCITY_BOUNDARY_CONDITIONS(Inclusion,
+                                                   "inclusion",
+                                                   "Implementation of the velocity boundary conditions for the "
+                                                   "``inclusion'' benchmark. See the manual and the Kronbichler, Heister "
+                                                   "and Bangerth paper on ASPECT for more information about this "
+                                                   "benchmark.")
+    }
+  }
 }
diff --git a/source/postprocess/duretz_et_al.cc b/source/postprocess/duretz_et_al.cc
index b802218..324dff5 100644
--- a/source/postprocess/duretz_et_al.cc
+++ b/source/postprocess/duretz_et_al.cc
@@ -2842,42 +2842,7 @@ namespace aspect
       }
 
 
-      // based on http://geodynamics.org/hg/cs/AMR/Discontinuous_Stokes with permission
-      void _Inclusion(double pos[2], double r_inclusion, double eta, double *vx, double *vy, double *p)
-      {
-        const double min_eta = 1.0;
-        const double max_eta = eta;
-        const double epsilon = 1; //strain rate
-        const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
-        std::complex<double> phi, psi, dphi;
-        const double offset[2]= {1.0, 1.0};
-        double r2_inclusion = r_inclusion * r_inclusion;
-
-        double x = pos[0]-offset[0];
-        double y = pos[1]-offset[1];
-        double r2 = x*x+y*y;
-
-        std::complex<double> z(x,y);
-        if (r2<r2_inclusion)
-          {
-            //inside the inclusion
-            phi=0;
-            dphi=0;
-            psi=-4*epsilon*(max_eta*min_eta/(min_eta+max_eta))*z;
-          }
-        else
-          {
-            //outside the inclusion
-            phi=-2*epsilon*A*r2_inclusion/z;
-            dphi=-phi/z;
-            psi=-2*epsilon*(min_eta*z+A*r2_inclusion*r2_inclusion/(z*z*z));
-          }
-        double visc = (r2<r2_inclusion)? max_eta : 1.0;
-        std::complex<double> v = (phi - z*conj(dphi) - conj(psi))/(2.0*visc);
-        *vx=v.real();
-        *vy=v.imag();
-        *p=-2*epsilon*dphi.real();
-      }
+
     }
   }
 }