Add variable wall props for LELAS and STRUCT 3D

Code/Source/svFSI/DISTRIBUTE.f

@@ -186,6 +186,10 @@ SUBROUTINE DISTRIBUTE
          CALL cm%bcast(cmmVarWall)
          CALL cm%bcast(shlEq)
          CALL cm%bcast(pstEq)
+!        Varwall properties---------------------------------------------
+         CALL cm%bcast(useVarWall)
+         CALL cm%bcast(nvwp)
+!        ---------------------------------------------------------------
          CALL cm%bcast(sstEq)
          CALL cm%bcast(cepEq)
          IF (rmsh%isReqd) THEN
@@ -256,6 +260,24 @@ SUBROUTINE DISTRIBUTE
          DEALLOCATE(tmpX)
       END IF
 
+!     Varwall properties------------------------------------------------
+!     Distribute variable wall properties (vWP0) to processors
+      flag = ALLOCATED(vWP0)
+      CALL cm%bcast(flag)
+      IF (flag) THEN
+         IF (cm%mas()) THEN
+            ALLOCATE(tmpX(nvwp,gtnNo))
+            tmpX = vWP0
+            DEALLOCATE(vWP0)
+         ELSE
+            ALLOCATE(tmpX(0,0))
+         END IF
+         ALLOCATE(vWP0(nvwp,tnNo))
+         vWP0 = LOCAL(tmpX)
+         DEALLOCATE(tmpX)
+      END IF
+!     ------------------------------------------------------------------
+
 !     Distribute initial flow quantities to processors
       flag = ALLOCATED(Pinit)
       CALL cm%bcast(flag)

Code/Source/svFSI/FSI.f

@@ -52,7 +52,7 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
       INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
       REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), yl(:,:),
      2   dl(:,:), bfl(:,:), fN(:,:), pS0l(:,:), pSl(:), ya_l(:),
-     3   lR(:,:), lK(:,:,:), lKd(:,:,:)
+     3   lR(:,:), lK(:,:,:), lKd(:,:,:), lVWP(:,:)
       REAL(KIND=RKIND), ALLOCATABLE :: xwl(:,:), xql(:,:), Nwx(:,:),
      2   Nwxx(:,:), Nqx(:,:)
 
@@ -72,7 +72,7 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
       ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN),
      2   dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN),
      3   pSl(nsymd), ya_l(eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN),
-     4   lKd(dof*nsd,eNoN,eNoN))
+     4   lKd(dof*nsd,eNoN,eNoN), lVWP(nvwp,eNoN))
 
 !     Loop over all elements of mesh
       DO e=1, lM%nEl
@@ -98,6 +98,12 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
             yl(:,a)  = Yg(:,Ac)
             dl(:,a)  = Dg(:,Ac)
             bfl(:,a) = Bf(:,Ac)
+
+!           Varwall properties -----------------------------------------
+!           Calculate local wall property
+            IF (useVarWall) lVWP(:,a) = vWP0(:,Ac)
+!           ------------------------------------------------------------
+
             IF (ALLOCATED(lM%fN)) THEN
                DO iFn=1, nFn
                   fN(:,iFn) = lM%fN((iFn-1)*nsd+1:iFn*nsd,e)
@@ -159,11 +165,11 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg)
 
                CASE (phys_lElas)
                   CALL LELAS3D(fs(1)%eNoN, w, fs(1)%N(:,g), Nwx, al, dl,
-     2               bfl, pS0l, pSl, lR, lK)
+     2               bfl, pS0l, pSl, lR, lK, lVWP)
 
                CASE (phys_struct)
                   CALL STRUCT3D(fs(1)%eNoN, nFn, w, fs(1)%N(:,g), Nwx,
-     2               al, yl, dl, bfl, fN, pS0l, pSl, ya_l, lR, lK)
+     2               al, yl, dl, bfl, fN, pS0l, pSl, ya_l, lR, lK, lVWP)
 
                CASE (phys_ustruct)
                   CALL USTRUCT3D_M(vmsStab, fs(1)%eNoN, fs(2)%eNoN, nFn,

Code/Source/svFSI/INITIALIZE.f

@@ -154,6 +154,7 @@ SUBROUTINE INITIALIZE(timeP)
       IF (dFlag) i = 3*tDof
       IF (pstEq) i = i + nsymd
       IF (sstEq) i = i + nsd
+!      IF (useVarWall) i = i + nvwp
       IF (cepEq) THEN
          i = i + nXion
          IF (cem%cpld) i = i + 1
@@ -663,6 +664,10 @@ SUBROUTINE FINALIZE
       IF (ALLOCATED(pSn))      DEALLOCATE(pSn)
       IF (ALLOCATED(pSa))      DEALLOCATE(pSa)
 
+!     Varwall properties -----------------------------------------------
+      IF (ALLOCATED(vWP0))      DEALLOCATE(vWP0)
+!     ------------------------------------------------------------------
+
       IF (ALLOCATED(Pinit))    DEALLOCATE(Pinit)
       IF (ALLOCATED(Vinit))    DEALLOCATE(Vinit)
       IF (ALLOCATED(Dinit))    DEALLOCATE(Dinit)

Code/Source/svFSI/LELAS.f

@@ -47,14 +47,17 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
 
       INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
       REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), dl(:,:),
-     2   bfl(:,:), pS0l(:,:), pSl(:), N(:), Nx(:,:), lR(:,:), lK(:,:,:)
+     2   bfl(:,:), pS0l(:,:), pSl(:), N(:), Nx(:,:), lR(:,:), lK(:,:,:),
+     3   lVWP(:,:)
 
       eNoN = lM%eNoN
 
 !     LELAS: dof = nsd
       ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), dl(tDof,eNoN),
      2   bfl(nsd,eNoN), pS0l(nsymd,eNoN), pSl(nsymd), N(eNoN),
-     3   Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN))
+     3   Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN),
+     4   lVWP(nvwp,eNoN))
+
 
 !     Loop over all elements of mesh
       DO e=1, lM%nEl
@@ -76,8 +79,15 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
             dl(:,a)  = Dg(:,Ac)
             bfl(:,a) = Bf(:,Ac)
             IF (ALLOCATED(pS0)) pS0l(:,a) = pS0(:,Ac)
+!           Varwall properties------------------------------------------
+!           Calculate local wall property
+            IF (useVarWall) THEN
+               lVWP(:,a) = vWP0(:,Ac)
+            END IF
+!        ---------------------------------------------------------------
          END DO
 
+
 !        Gauss integration
          lR = 0._RKIND
          lK = 0._RKIND
@@ -92,7 +102,7 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
             pSl = 0._RKIND
             IF (nsd .EQ. 3) THEN
                CALL LELAS3D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, lR,
-     2            lK)
+     2            lK, lVWP)
 
             ELSE IF (nsd .EQ. 2) THEN
                CALL LELAS2D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, lR,
@@ -128,39 +138,33 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg)
       END SUBROUTINE CONSTRUCT_LELAS
 !####################################################################
       PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
-     2   lR, lK)
+     2   lR, lK, lVWP)
       USE COMMOD
       IMPLICIT NONE
       INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
       REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN),
-     2   al(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN), pS0l(6,eNoN)
+     2   al(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN), pS0l(6,eNoN),
+     3   lVWP(nvwp,eNoN)
       REAL(KIND=RKIND), INTENT(INOUT) :: pSl(6), lR(dof,eNoN),
      2   lK(dof*dof,eNoN,eNoN)
 
       INTEGER(KIND=IKIND) a, b, i, j, k
       REAL(KIND=RKIND) NxdNx, rho, elM, nu, lambda, mu, divD, T1, amd,
-     2   wl, lDm, ed(6), ud(3), f(3), S0(6), S(6)
+     2   wl, lDm, ed(6), ud(3), f(3), S0(6), S(6), eVWP(nvwp), Cst(6,6)
 
       rho  = eq(cEq)%dmn(cDmn)%prop(solid_density)
-      elM  = eq(cEq)%dmn(cDmn)%prop(elasticity_modulus)
-      nu   = eq(cEq)%dmn(cDmn)%prop(poisson_ratio)
       f(1) = eq(cEq)%dmn(cDmn)%prop(f_x)
       f(2) = eq(cEq)%dmn(cDmn)%prop(f_y)
       f(3) = eq(cEq)%dmn(cDmn)%prop(f_z)
       i    = eq(cEq)%s
       j    = i + 1
       k    = j + 1
 
-      lambda = elM*nu / (1._RKIND+nu) / (1._RKIND-2._RKIND*nu)
-      mu     = elM * 0.5_RKIND / (1._RKIND+nu)
-      lDm    = lambda/mu
-      T1     = eq(cEq)%af*eq(cEq)%beta*dt*dt
-      amd    = eq(cEq)%am/T1*rho
-      wl     = w*T1*mu
-
       ed = 0._RKIND
       ud = -f
       S0 = 0._RKIND
+      eVWP = 0._RKIND
+
       DO a=1, eNoN
          ud(1) = ud(1) + N(a)*(al(i,a)-bfl(1,a))
          ud(2) = ud(2) + N(a)*(al(j,a)-bfl(2,a))
@@ -173,22 +177,49 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
          ed(5) = ed(5) + Nx(3,a)*dl(j,a) + Nx(2,a)*dl(k,a)
          ed(6) = ed(6) + Nx(1,a)*dl(k,a) + Nx(3,a)*dl(i,a)
 
+!     ------------------------------------------------------------------
+!     Calculate local wall property
+!     Don't use if calculating mesh
+         IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN
+            eVWP(:) = eVWP(:) + N(a)*lVWP(:,a)
+         END IF
+!     ------------------------------------------------------------------
+
          S0(1) = S0(1) + N(a)*pS0l(1,a)
          S0(2) = S0(2) + N(a)*pS0l(2,a)
          S0(3) = S0(3) + N(a)*pS0l(3,a)
          S0(4) = S0(4) + N(a)*pS0l(4,a)
          S0(5) = S0(5) + N(a)*pS0l(5,a)
          S0(6) = S0(6) + N(a)*pS0l(6,a)
       END DO
+
+
+      IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN
+         elM  = eVWP(1)
+         nu   = evWP(2)
+      ELSE
+         elM  = eq(cEq)%dmn(cDmn)%prop(elasticity_modulus)
+         nu   = eq(cEq)%dmn(cDmn)%prop(poisson_ratio)
+      END IF
+
+      lambda = elM*nu / (1._RKIND+nu) / (1._RKIND-2._RKIND*nu)
+      mu     = elM * 0.5_RKIND / (1._RKIND+nu)
+      lDm    = lambda/mu
+      T1     = eq(cEq)%af*eq(cEq)%beta*dt*dt
+      amd    = eq(cEq)%am/T1*rho
+      wl     = w*T1*mu
+
       divD = lambda*(ed(1) + ed(2) + ed(3))
 
+
 !     Stress in Voigt notation
       S(1) = divD + 2._RKIND*mu*ed(1)
       S(2) = divD + 2._RKIND*mu*ed(2)
       S(3) = divD + 2._RKIND*mu*ed(3)
       S(4) = mu*ed(4)  ! 2*eps_12
       S(5) = mu*ed(5)  ! 2*eps_23
       S(6) = mu*ed(6)  ! 2*eps_13
+
       pSl  = S
 
 !     Add prestress contribution
@@ -253,7 +284,7 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
 
       INTEGER(KIND=IKIND) a, b, i, j
       REAL(KIND=RKIND) NxdNx, rho, elM, nu, lambda, mu, divD, T1, amd,
-     2   wl, lDm, ed(3), ud(2), f(2), S0(3), S(3)
+     2   T2, wl, lDm, ed(3), ud(2), f(2), S0(3), S(3)
 
       rho  = eq(cEq)%dmn(cDmn)%prop(solid_density)
       elM  = eq(cEq)%dmn(cDmn)%prop(elasticity_modulus)
@@ -296,6 +327,8 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
 !     Add prestress contribution
       S = pSl + S0
 
+!     Need to add variable wall tangent matrix
+
       DO a=1, eNoN
          lR(1,a) = lR(1,a) + w*(rho*N(a)*ud(1) + Nx(1,a)*S(1) +
      2      Nx(2,a)*S(3))
@@ -306,9 +339,9 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
          DO b=1, eNoN
             NxdNx = Nx(1,a)*Nx(1,b) + Nx(2,a)*Nx(2,b)
 
-            T1 = amd*N(a)*N(b)/mu + NxdNx
+            T2 = amd*N(a)*N(b)/mu + NxdNx
 
-            lK(1,a,b) = lK(1,a,b) + wl*(T1
+            lK(1,a,b) = lK(1,a,b) + wl*(T2
      2         + (1._RKIND + lDm)*Nx(1,a)*Nx(1,b))
 
             lK(2,a,b) = lK(2,a,b) + wl*(lDm*Nx(1,a)*Nx(2,b)
@@ -317,7 +350,7 @@ PURE SUBROUTINE LELAS2D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl,
             lK(dof+1,a,b) = lK(dof+1,a,b) + wl*(lDm*Nx(2,a)*Nx(1,b)
      2         + Nx(1,a)*Nx(2,b))
 
-            lK(dof+2,a,b) = lK(dof+2,a,b) + wl*(T1
+            lK(dof+2,a,b) = lK(dof+2,a,b) + wl*(T2
      2         + (1._RKIND + lDm)*Nx(2,a)*Nx(2,b))
          END DO
       END DO

Code/Source/svFSI/MATMODELS.f

@@ -38,14 +38,15 @@
 
 !     Compute 2nd Piola-Kirchhoff stress and material stiffness tensors
 !     including both dilational and isochoric components
-      SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm)
+      SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm, eVWP)
       USE MATFUN
       USE COMMOD
       IMPLICIT NONE
       TYPE(dmnType), INTENT(IN) :: lDmn
       INTEGER(KIND=IKIND), INTENT(IN) :: nfd
       REAL(KIND=RKIND), INTENT(IN) :: F(nsd,nsd), fl(nsd,nfd), ya
       REAL(KIND=RKIND), INTENT(OUT) :: S(nsd,nsd), Dm(nsymd,nsymd)
+      REAL(KIND=RKIND), INTENT(IN), OPTIONAL :: eVWP(nvwp)
 
       TYPE(stModelType) :: stM
       REAL(KIND=RKIND) :: nd, Kp, J, J2d, J4d, trE, p, pl, Inv1, Inv2,
@@ -129,7 +130,12 @@ SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm)
 
 !     NeoHookean model
       CASE (stIso_nHook)
-         g1 = 2._RKIND * stM%C10
+         IF (useVarWall) THEN
+!           Converting elastic modulus and poisson ratio to g1
+            g1 = eVWP(1)* 0.5_RKIND/(1._RKIND+eVWP(2))
+         ELSE
+            g1 = 2._RKIND * stM%C10
+         END IF
          Sb = g1*IDm
 
 !        Fiber reinforcement/active stress

Code/Source/svFSI/MOD.f

@@ -837,6 +837,8 @@ MODULE COMMOD
       LOGICAL cmmInit
 !     Whether variable wall properties are used for CMM
       LOGICAL cmmVarWall
+!     Whether variable wall properties should be read and used
+      LOGICAL useVarWall
 !     Whether shell equation is being solved
       LOGICAL shlEq
 !     Whether PRESTRESS is being solved
@@ -893,6 +895,8 @@ MODULE COMMOD
       INTEGER(KIND=IKIND) rsTS
 !     Number of stress values to be stored
       INTEGER(KIND=IKIND) nsymd
+!     Number of variable wall properties to read in from mesh
+      INTEGER(KIND=IKIND) nvwp
 
 !     REAL VARIABLES
 !     Time step size
@@ -933,7 +937,7 @@ MODULE COMMOD
       REAL(KIND=RKIND), ALLOCATABLE :: Ao(:,:)
 !     New time derivative of variables
       REAL(KIND=RKIND), ALLOCATABLE :: An(:,:)
-!     Old integrated variables (dissplacement)
+!     Old integrated variables (displacement)
       REAL(KIND=RKIND), ALLOCATABLE :: Do(:,:)
 !     New integrated variables
       REAL(KIND=RKIND), ALLOCATABLE :: Dn(:,:)
@@ -964,6 +968,9 @@ MODULE COMMOD
       REAL(KIND=RKIND), ALLOCATABLE :: pSn(:,:)
       REAL(KIND=RKIND), ALLOCATABLE :: pSa(:)
 
+!     Variables for variable wall properties
+      REAL(KIND=RKIND), ALLOCATABLE :: vWP0(:,:)
+
 !     Temporary storage for initializing state variables
       REAL(KIND=RKIND), ALLOCATABLE :: Pinit(:)
       REAL(KIND=RKIND), ALLOCATABLE :: Vinit(:,:)