From cfd10ed9612c7a0aeb016d62e32e9dd4380ce4c1 Mon Sep 17 00:00:00 2001
From: Matthias Schmiddunser <matthias.schmiddunser@uni-tuebingen.de>
Date: Fri, 19 Sep 2014 18:14:20 +0200
Subject: [PATCH] Addition temperature gradient at bottom

5th temperature option added:
Radial Gradient from (0,y0,0) to indenter edge
Updated input template
---
 src/define_bc_rot_subduction.f90 |  9 ++++-
 test/input_schmiddu.txt          | 58 +++++++++++++++++---------------
 2 files changed, 38 insertions(+), 29 deletions(-)

diff --git a/src/define_bc_rot_subduction.f90 b/src/define_bc_rot_subduction.f90
index 1db5bc59..956f880f 100644
--- a/src/define_bc_rot_subduction.f90
+++ b/src/define_bc_rot_subduction.f90
@@ -609,7 +609,7 @@ do i=1,osolve%nnode
             select case(tempcase)
                 case (1)
                     osolve%temp(i) = temp2+osolve%x(i)*(1.0d0-temp2)
-                case (2,3,4)
+                case (2,3,4,5)
                     osolve%temp(i)=1.d0
                 case default
                     osolve%temp(i)=1.d0
@@ -631,6 +631,8 @@ do i=1,osolve%nnode
                     osolve%temp(i) = 1.d0+osolve%x(i)/rx*(temp2-1.0d0)
                 case (4)
                     osolve%temp(i) = 1.d0+osolve%x(i)/yfix_xedge*(temp2-1.0d0)
+                case (5)
+                    osolve%temp(i) = 1.d0+distz0plane/rx*(temp2-1.0d0)
                 case default
                     osolve%temp(i)=1.d0                            
             end select
@@ -663,6 +665,11 @@ do i=1,osolve%nnode
                     tempstart = 1.d0+xminedge/yfix_xedge*(temp2-1.0d0)
                     tempend = 1.d0
                     osolve%temp(i)=tempstart+distzonerad*(tempend-tempstart)
+                case (5)
+                    distz0plane = sqrt(xedge**2+((yedge-y0)*rx/(0.5d0*wy))**2)
+                    tempstart = 1.d0+distz0plane/rx*(temp2-1.0d0)
+                    tempend = 1.d0
+                    osolve%temp(i)=tempstart+distzonerad*(tempend-tempstart)
                 case default
                     osolve%temp(i)=1.d0                            
             end select 
diff --git a/test/input_schmiddu.txt b/test/input_schmiddu.txt
index a7c70aa7..cfae8011 100644
--- a/test/input_schmiddu.txt
+++ b/test/input_schmiddu.txt
@@ -428,34 +428,36 @@ Types requiring additional input values are listed below.
                       subtract 1.0 from all u component velocities along all
                       faces.
 
-- rot_subduction   -> upper section of a downwards rotating squashed sphere or cylinder
-                          (rotation around the -y-axis)
-                   -> [dp]bc_param1 is y0 the center of the rotating element
-                   -> [dp]bc_param2 is h, the maximum height of the rotating element
-                          influx within element will be mass-balanced by bottom outflux
-                   -> [dp]bc_param3 is rx, the x-width of the rotating element
-                   -> [dp]bc_param4 is wy, the y-width of the rotating element. wy>0 will 
-                          produce a squashed sphere, wy<0 a cylinder 
-                   -> [dp]bc_param5 is vin, velocity of the rotating element at (0,y0,h)
-                   -> [dp]bc_param6 is vtop, the velocity of influx above the rotating element
-                          (ususally vin)
-                   -> [dp]bc_param7 is vback, the velocity at the rights side and along 
-                          the other parts of bottom and left side
-                          NOTE: measured in -y direction
-                   -> [dp]bc_param8 is nelemx the number of elements across
-                          which to spread the velocity BCs in x- or y-direction 
-                   -> [dp]bc_param9 is nelemz the number of elements across
-                          which to spread the velocity BCs on the z-direction
-                   -> [dp]bc_param10 is tempcase, the selection for temperature boundary
-                          conditions at the base
-                          0: uniform temperature
-                          1: uniform temp2 for indenter footprint, tempscale elsewhere
-                          2: overall gradient from temp2 at x=0 to tepmscale at x=1
-                          3: gradient from tempscale at x=0 to temp2 at x=rx within the 
-                             indenter footprint, tempscale elsewhere
-                          4: gradient from tempscale at x=0 to temp2 at the edge of the 
-                             rotating element, tempscale elsewhere
-                   -> [dp]bc_param11 is temp2 in Celsius and should be less than tempscale                  
+- rot_subduction   
+    -> upper section of a downwards rotating squashed sphere or cylinder 
+        (rotation around the -y-axis) 
+    -> [dp]bc_param1 is y0 the center of the rotating element 
+    -> [dp]bc_param2 is h, the maximum height of the rotating element influx 
+        within element will be mass-balanced by bottom outflux 
+    -> [dp]bc_param3 is rx, the x-width of the rotating element 
+    -> [dp]bc_param4 is wy, the y-width of the rotating element.  wy>0 will 
+        produce a squashed sphere, wy<0 a cylinder 
+    -> [dp]bc_param5 is vin, velocity of the rotating element at (0,y0,h) 
+    -> [dp]bc_param6 is vtop, the velocity of influx above the rotating 
+        element (ususally vin) -> [dp]bc_param7 is vback, the velocity at the 
+        rights side and along the other parts of bottom and left side
+        NOTE: measured in -y direction 
+    -> [dp]bc_param8 is nelemx the number of elements across which to spread 
+        the velocity BCs in x- or y-direction 
+    -> [dp]bc_param9 is nelemz the number of elements across which to spread 
+        the velocity BCs on the z-direction 
+    -> [dp]bc_param10 is tempcase, the selection for temperature boundary 
+        conditions at the base 
+            0: uniform temperature 
+            1: uniform temp2 for indenter footprint, tempscale elsewhere
+            2: overall gradient from temp2 at x=0 to tepmscale at x=1
+            3: gradient from tempscale at x=0 to temp2 at x=rx within the 
+               indenter footprint, tempscale elsewhere
+            4: gradient from tempscale at x=0 to temp2 at the edge of the 
+               rotating element, tempscale elsewhere
+            5: radial gradient from tempscale at (x=0,y=y0) to temp 2 at the 
+               edge of the indenter, tempscale elsewhere
+    -> [dp]bc_param11 is temp2 in Celsius and should be less than tempscale                  
                       
 - iso_only         -> Applies zero displacement boundary conditions to all faces
                       and dynamically calculates an initial optimal time step
-- 
GitLab