[cig-commits] [commit] master: move more inclusion stuff out (f789e0a)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed May 21 04:23:46 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/c883fc069e0dcc93cce45e1d52c84c05c9c4385b...1f4ab820b0191d9a4058b78c91b629bc94c437ec
>---------------------------------------------------------------
commit f789e0a380d8763a55842ce0938ef3d9e572a926
Author: Timo Heister <timo.heister at gmail.com>
Date: Tue May 20 11:38:26 2014 -0400
move more inclusion stuff out
>---------------------------------------------------------------
f789e0a380d8763a55842ce0938ef3d9e572a926
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();
- }
+
}
}
}
More information about the CIG-COMMITS
mailing list