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(:,:)