adding pre_stress hardcode patch

Author: mcarcana

src/pre_process/include/3dHardcodedIC.fpp

@@ -4,6 +4,11 @@
     real(kind(0d0)) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph
 
     real(kind(0d0)) :: eps
+   
+    real(kind(0d0)) :: rcoord, theta, phi, xi_sph, x_bcen, y_bcen, z_bcen, Rinit
+    real(kind(0d0)) :: x_ccs, y_ccs, z_ccs
+    real(kind(0d0)), dimension(num_dims) :: xi_cart
+    integer :: l
 
     eps = 1e-9
 #:enddef
@@ -56,6 +61,33 @@
             q_prim_vf(advxe)%sf(i, j, k) = patch_icpp(1)%alpha(2)
         end if
 
+    case (302) ! pre_stress for hyperelasticity, bubble in material
+        R0ref = 30E-6    ! equilibrium radius
+        Rinit = patch_icpp(3)%radius ! initial radius
+        x_bcen = patch_icpp(3)%x_centroid
+        y_bcen = patch_icpp(3)%y_centroid
+        z_bcen = patch_icpp(3)%z_centroid
+        x_ccs = x_cc(i) - x_bcen
+        y_ccs = y_cc(j) - y_bcen
+        z_ccs = z_cc(k) - z_bcen
+        rcoord = sqrt(x_ccs**2 + y_ccs**2 + z_ccs**2)
+        phi = atan2(y_ccs, x_ccs)
+        theta = atan2(sqrt(x_ccs**2 + y_ccs**2), z_ccs)
+        !spherical coord, assuming Rmax=1
+        xi_sph = (rcoord**3 - R0ref**3 + Rinit**3)**(1d0/3d0) 
+        xi_cart(1) = xi_sph*sin(theta)*cos(phi)
+        xi_cart(2) = xi_sph*sin(theta)*sin(phi)
+        xi_cart(3) = xi_sph*cos(theta)
+        ! shift back
+        xi_cart(1) = xi_cart(1) + x_bcen
+        xi_cart(2) = xi_cart(2) + y_bcen
+        xi_cart(3) = xi_cart(3) + z_bcen
+!        print *, 'xi_cart(1) ::', xi_cart(1), 'xi_cart(2) ::', xi_cart(2), 'xi_cart(3) ::', xi_cart(3)
+        ! assigning the reference map to the q_prim vector field
+        do l = 1, 3
+            q_prim_vf(l + xibeg - 1)%sf(i, j, k) = xi_cart(l)
+        end do
+
         ! Put your variable assignments here
     case default
         call s_int_to_str(patch_id, iStr)

src/pre_process/m_assign_variables.fpp

@@ -500,29 +500,20 @@ contains
         end if
 
         ! Elastic Shear Stress
-        if (hyperelasticity) then
-
-            if (pre_stress) then ! pre stressed initial condition in spatial domain
-                rcoord = sqrt((x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2))
-                theta = atan2(y_cc(k), x_cc(j))
-                phi = atan2(sqrt(x_cc(j)**2 + y_cc(k)**2), z_cc(l))
-                !spherical coord, assuming Rmax=1
-                xi_sph = (rcoord**3 - R0ref**3 + 1d0)**(1d0/3d0)
-                xi_cart(1) = xi_sph*sin(phi)*cos(theta)
-                xi_cart(2) = xi_sph*sin(phi)*sin(theta)
-                xi_cart(3) = xi_sph*cos(phi)
-            else
-                xi_cart(1) = x_cc(j)
-                xi_cart(2) = y_cc(k)
-                xi_cart(3) = z_cc(l)
+        if (hyperelasticity .and. .not. pre_stress) then
+            xi_cart(1) = x_cc(j)
+            if (p > 0) then
+               xi_cart(2) = y_cc(k)
+               xi_cart(3) = z_cc(l)
+            elseif (n > 0) then
+               xi_cart(2) = y_cc(k)
             end if
 
             ! assigning the reference map to the q_prim vector field
             do i = 1, num_dims
-                q_prim_vf(i + xibeg - 1)%sf(j, k, l) = eta*xi_cart(i) + &
+               q_prim_vf(i + xibeg - 1)%sf(j, k, l) = eta*xi_cart(i) + &
                                                        (1d0 - eta)*orig_prim_vf(i + xibeg - 1)
             end do
-
         end if
 
         if (mpp_lim .and. bubbles) then