From 533f8c90c929a6e420db664de77c70568b9e196c Mon Sep 17 00:00:00 2001 From: Erica Date: Fri, 23 Jul 2021 16:25:42 -0700 Subject: [PATCH 01/10] Single processor varwall nHk --- Code/Source/svFSI/MATMODELS.f | 16 ++++++++++++++-- Code/Source/svFSI/MOD.f | 12 ++++++++++++ Code/Source/svFSI/POST.f | 24 ++++++++++++++++++++---- Code/Source/svFSI/READMSH.f | 34 +++++++++++++++++++++++++++++++++- Code/Source/svFSI/STRUCT.f | 32 ++++++++++++++++++++++++-------- Code/Source/svFSI/VTKXML.f | 2 +- 6 files changed, 104 insertions(+), 16 deletions(-) diff --git a/Code/Source/svFSI/MATMODELS.f b/Code/Source/svFSI/MATMODELS.f index 36733867..0eaf7497 100755 --- a/Code/Source/svFSI/MATMODELS.f +++ b/Code/Source/svFSI/MATMODELS.f @@ -38,7 +38,7 @@ ! 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 @@ -46,6 +46,9 @@ SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm) 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) +! VARIABLE WALL PROPERTIES - SCHWARZ JULY 2021 --------------------- + 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 +132,16 @@ 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)) +! PRINT *, "USING VARWALL" +! PRINT *, g1 + ELSE + g1 = 2._RKIND * stM%C10 +! PRINT *, "NOT USING VARWALL" +! PRINT *, g1 + END IF Sb = g1*IDm ! Fiber reinforcement/active stress diff --git a/Code/Source/svFSI/MOD.f b/Code/Source/svFSI/MOD.f index 743f25ad..ebdc78be 100755 --- a/Code/Source/svFSI/MOD.f +++ b/Code/Source/svFSI/MOD.f @@ -837,6 +837,10 @@ MODULE COMMOD LOGICAL cmmInit ! Whether variable wall properties are used for CMM LOGICAL cmmVarWall +! SCHWARZ July 2021---------------------------------------------- +! 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 +897,10 @@ MODULE COMMOD INTEGER(KIND=IKIND) rsTS ! Number of stress values to be stored INTEGER(KIND=IKIND) nsymd +! SCHWARZ July 2021---------------------------------------------- +! Number of variable wall properties to read in from mesh + INTEGER(KIND=IKIND) nvwp +! --------------------------------------------------------------- ! REAL VARIABLES ! Time step size @@ -964,6 +972,10 @@ MODULE COMMOD REAL(KIND=RKIND), ALLOCATABLE :: pSn(:,:) REAL(KIND=RKIND), ALLOCATABLE :: pSa(:) +! Variables for variable wall properties - SCHWARZ July 2021---- + REAL(KIND=RKIND), ALLOCATABLE :: vWP0(:,:) +! -------------------------------------------------------------- + ! Temporary storage for initializing state variables REAL(KIND=RKIND), ALLOCATABLE :: Pinit(:) REAL(KIND=RKIND), ALLOCATABLE :: Vinit(:,:) diff --git a/Code/Source/svFSI/POST.f b/Code/Source/svFSI/POST.f index cf1c7568..e89d45e9 100755 --- a/Code/Source/svFSI/POST.f +++ b/Code/Source/svFSI/POST.f @@ -521,12 +521,13 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) REAL(KIND=RKIND) w, Jac, detF, Je, ya, Ja, elM, nu, lambda, mu, 2 p, trS, vmises, xi(nsd), xi0(nsd), xp(nsd), ed(nsymd), 3 Im(nsd,nsd), F(nsd,nsd), C(nsd,nsd), Eg(nsd,nsd), P1(nsd,nsd), - 4 S(nsd,nsd), sigma(nsd,nsd), Dm(nsymd,nsymd) + 4 S(nsd,nsd), sigma(nsd,nsd), Dm(nsymd,nsymd), eVWP(nvwp) TYPE(fsType) :: fs INTEGER, ALLOCATABLE :: eNds(:) REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), dl(:,:), yl(:,:), - 2 fN(:,:), resl(:), Nx(:,:), N(:), sA(:), sF(:,:), sE(:) + 2 fN(:,:), resl(:), Nx(:,:), N(:), sA(:), sF(:,:), + 3 sE(:), lVWP(:,:) dof = eq(iEq)%dof i = eq(iEq)%s @@ -555,7 +556,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) ALLOCATE (sA(tnNo), sF(m,tnNo), sE(lM%nEl), xl(nsd,fs%eNoN), 2 dl(tDof,fs%eNoN), yl(tDof,fs%eNoN), fN(nsd,nFn), resl(m), - 3 Nx(nsd,fs%eNoN), N(fs%eNoN)) + 3 Nx(nsd,fs%eNoN), N(fs%eNoN), lVWP(nvwp,fs%eNoN)) sA = 0._RKIND sF = 0._RKIND @@ -594,6 +595,12 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) xl(:,a) = x(:,Ac) dl(:,a) = lD(:,Ac) yl(:,a) = lY(:,Ac) +! Variable wall - SCHWARZ July 2021--------------------------- +! Calculate local wall property + IF (useVarWall) THEN + lVWP(:,a) = vWP0(:,Ac) + END IF +! ------------------------------------------------------------ END DO Je = 0._RKIND @@ -607,6 +614,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) Im = MAT_ID(nsd) F = Im + eVWP = 0._RKIND DO a=1, fs%eNoN IF (nsd .EQ. 3) THEN F(1,1) = F(1,1) + Nx(1,a)*dl(i,a) @@ -618,6 +626,13 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) F(3,1) = F(3,1) + Nx(1,a)*dl(k,a) F(3,2) = F(3,2) + Nx(2,a)*dl(k,a) F(3,3) = F(3,3) + Nx(3,a)*dl(k,a) +! Variable wall - SCHWARZ July 2021--------------------- +! Calculate local wall property + IF (useVarWall) THEN + eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) + END IF +! ------------------------------------------------------ + ELSE F(1,1) = F(1,1) + Nx(1,a)*dl(i,a) F(1,2) = F(1,2) + Nx(2,a)*dl(i,a) @@ -734,7 +749,8 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) IF (.NOT.ISZERO(detF)) sigma(:,:) = sigma(:,:) / detF ELSE IF (cPhys .EQ. phys_struct) THEN - CALL GETPK2CC(eq(iEq)%dmn(cDmn), F, nFn, fN, ya, S,Dm) + CALL GETPK2CC(eq(iEq)%dmn(cDmn), F, nFn, fN, ya, + 2 S,Dm, eVWP) P1 = MATMUL(F, S) sigma = MATMUL(P1, TRANSPOSE(F)) IF (.NOT.ISZERO(detF)) sigma(:,:) = sigma(:,:) / detF diff --git a/Code/Source/svFSI/READMSH.f b/Code/Source/svFSI/READMSH.f index a56269a6..c0177cc2 100755 --- a/Code/Source/svFSI/READMSH.f +++ b/Code/Source/svFSI/READMSH.f @@ -50,7 +50,7 @@ SUBROUTINE READMSH(list) INTEGER(KIND=IKIND) :: i, j, iM, iFa, a, b, Ac, e, lDof, lnNo REAL(KIND=RKIND) :: maxX(nsd), minX(nsd), fibN(nsd), rtmp CHARACTER(LEN=stdL) :: ctmp, fExt - TYPE(listType), POINTER :: lPtr, lPM + TYPE(listType), POINTER :: lPtr, lPM, lPtr2 TYPE(stackType) :: avNds TYPE(fileType) :: fTmp @@ -408,6 +408,38 @@ SUBROUTINE READMSH(list) END DO END IF +! Read variable wall properties - SCHWARZ July 2021 ---------------- + flag = .FALSE. + DO iM=1, nMsh + lPM => list%get(msh(iM)%name,"Add mesh",iM) + lPtr2 => lPM%get(nvwp,"Number of variable wall properties") + lPtr => lPM%get(cTmp, "Variable wall properties file path") + IF (ASSOCIATED(lPtr) .AND. ASSOCIATED(lPtr2)) THEN + IF (rmsh%isReqd) THEN + err = "Variable wall properties "// + 2 "is not currently allowed with remeshing" + END IF + flag = .TRUE. + useVarWall = .TRUE. + ALLOCATE(msh(iM)%x(nvwp,msh(iM)%gnNo)) + msh(iM)%x = 0._RKIND + CALL READVTUPDATA(msh(iM), cTmp, "varWallProps", nvwp, 1) + END IF + END DO + IF (flag) THEN + ALLOCATE(vWP0(nvwp,gtnNo)) + vWP0 = 0._RKIND + DO iM=1, nMsh + IF (.NOT.ALLOCATED(msh(iM)%x)) CYCLE + DO a=1, msh(iM)%gnNo + Ac = msh(iM)%gN(a) + vWP0(:,Ac) = msh(iM)%x(:,a) + END DO + DEALLOCATE(msh(iM)%x) + END DO + END IF +! ------------------------------------------------------------------ + ! Load any initial data (velocity, pressure, displacements) IF (.NOT.resetSim) CALL LOADVARINI(list) diff --git a/Code/Source/svFSI/STRUCT.f b/Code/Source/svFSI/STRUCT.f index 933abcf5..0f6d7ac5 100755 --- a/Code/Source/svFSI/STRUCT.f +++ b/Code/Source/svFSI/STRUCT.f @@ -50,7 +50,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:) REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), yl(:,:), 2 dl(:,:), bfl(:,:), fN(:,:), pS0l(:,:), pSl(:), ya_l(:), N(:), - 3 Nx(:,:), lR(:,:), lK(:,:,:) + 3 Nx(:,:), lR(:,:), lK(:,:,:), lVWP(:,:) eNoN = lM%eNoN nFn = lM%nFn @@ -60,7 +60,7 @@ SUBROUTINE CONSTRUCT_dSOLID(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), N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), - 4 lK(dof*dof,eNoN,eNoN)) + 4 lK(dof*dof,eNoN,eNoN), lVWP(nvwp,eNoN)) ! Loop over all elements of mesh DO e=1, lM%nEl @@ -84,6 +84,13 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) yl(:,a) = Yg(:,Ac) dl(:,a) = Dg(:,Ac) bfl(:,a) = Bf(:,Ac) + +! Variable wall - SCHWARZ July 2021--------------------------- +! Calculate local wall property + IF (useVarWall) THEN + lVWP(:,a) = vWP0(:,Ac) + END IF +! ------------------------------------------------------------ IF (ALLOCATED(lM%fN)) THEN DO iFn=1, nFn fN(:,iFn) = lM%fN((iFn-1)*nsd+1:iFn*nsd,e) @@ -107,7 +114,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) pSl = 0._RKIND IF (nsd .EQ. 3) THEN CALL STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - 2 pS0l, pSl, ya_l, lR, lK) + 2 pS0l, pSl, ya_l, lR, lK, lVWP) ELSE IF (nsd .EQ. 2) THEN CALL STRUCT2D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, @@ -138,20 +145,20 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) END DO ! e: loop DEALLOCATE(ptr, xl, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, N, Nx, - 2 lR, lK) + 2 lR, lK, lVWP) RETURN END SUBROUTINE CONSTRUCT_dSOLID !#################################################################### SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - 2 pS0l, pSl, ya_l, lR, lK) + 2 pS0l, pSl, ya_l, lR, lK, lVWP) USE COMMOD USE ALLFUN IMPLICIT NONE INTEGER(KIND=IKIND), INTENT(IN) :: eNoN, nFn REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN), 2 al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN), - 3 fN(3,nFn), pS0l(6,eNoN), ya_l(eNoN) + 3 fN(3,nFn), pS0l(6,eNoN), ya_l(eNoN), lVWP(nvwp,eNoN) REAL(KIND=RKIND), INTENT(OUT) :: pSl(6) REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN), 2 lK(dof*dof,eNoN,eNoN) @@ -159,7 +166,7 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, INTEGER(KIND=IKIND) :: a, b, i, j, k REAL(KIND=RKIND) :: rho, dmp, T1, amd, afl, ya_g, fb(3), ud(3), 2 NxSNx, BmDBm, F(3,3), S(3,3), P(3,3), Dm(6,6), DBm(6,3), - 3 Bm(6,3,eNoN), S0(3,3) + 3 Bm(6,3,eNoN), S0(3,3), eVWP(nvwp) ! Define parameters rho = eq(cEq)%dmn(cDmn)%prop(solid_density) @@ -181,6 +188,8 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, F(3,3) = 1._RKIND S0 = 0._RKIND ya_g = 0._RKIND + eVWP = 0._RKIND + DO a=1, eNoN ud(1) = ud(1) + N(a)*(rho*(al(i,a)-bfl(1,a)) + dmp*yl(i,a)) ud(2) = ud(2) + N(a)*(rho*(al(j,a)-bfl(2,a)) + dmp*yl(j,a)) @@ -196,6 +205,13 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, F(3,2) = F(3,2) + Nx(2,a)*dl(k,a) F(3,3) = F(3,3) + Nx(3,a)*dl(k,a) +! Variable wall - SCHWARZ July 2021--------------------------------- +! Calculate local wall property + IF (useVarWall) THEN + eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) + END IF +! ------------------------------------------------------------------ + S0(1,1) = S0(1,1) + N(a)*pS0l(1,a) S0(2,2) = S0(2,2) + N(a)*pS0l(2,a) S0(3,3) = S0(3,3) + N(a)*pS0l(3,a) @@ -211,7 +227,7 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, ! 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in ! Voigt notationa (Dm) - CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm) + CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm, eVWP) ! Prestress pSl(1) = S(1,1) diff --git a/Code/Source/svFSI/VTKXML.f b/Code/Source/svFSI/VTKXML.f index 30317f43..7fbb80b0 100755 --- a/Code/Source/svFSI/VTKXML.f +++ b/Code/Source/svFSI/VTKXML.f @@ -255,7 +255,7 @@ SUBROUTINE READVTUPDATA(lM, fName, kwrd, m, idx) REAL(KIND=RKIND), ALLOCATABLE :: tmpR(:,:) iStat = 0 - std = " Loading file <"//TRIM(fName)//">" + std = " Loading file <"//TRIM(fName)//"> "//m CALL loadVTK(vtu, fName, iStat) IF (iStat .LT. 0) err = "VTU file read error (init)" From 1cf4d88f080db564857bb9f72ebd8fa50c654798 Mon Sep 17 00:00:00 2001 From: Erica Date: Mon, 26 Jul 2021 12:21:14 -0700 Subject: [PATCH 02/10] Adding multi-processor, error in DISTRIBUTE.f, though --- Code/Source/svFSI/DISTRIBUTE.f | 21 +++++++++++++++++++++ Code/Source/svFSI/INITIALIZE.f | 15 +++++++++++++++ Code/Source/svFSI/READFILES.f | 3 +++ Code/Source/svFSI/REMESH.f | 4 ++++ 4 files changed, 43 insertions(+) diff --git a/Code/Source/svFSI/DISTRIBUTE.f b/Code/Source/svFSI/DISTRIBUTE.f index 756c9c72..a13c4e20 100755 --- a/Code/Source/svFSI/DISTRIBUTE.f +++ b/Code/Source/svFSI/DISTRIBUTE.f @@ -186,6 +186,9 @@ SUBROUTINE DISTRIBUTE CALL cm%bcast(cmmVarWall) CALL cm%bcast(shlEq) CALL cm%bcast(pstEq) +! Variable wall properties - SCHWARZ July 2021------------------- + CALL cm%bcast(useVarWall) +! --------------------------------------------------------------- CALL cm%bcast(sstEq) CALL cm%bcast(cepEq) IF (rmsh%isReqd) THEN @@ -256,6 +259,24 @@ SUBROUTINE DISTRIBUTE DEALLOCATE(tmpX) END IF +! Variable wall properties - SCHWARZ July 2021 --------------------- +! 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) diff --git a/Code/Source/svFSI/INITIALIZE.f b/Code/Source/svFSI/INITIALIZE.f index 7cfe05c2..c35fb417 100755 --- a/Code/Source/svFSI/INITIALIZE.f +++ b/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 @@ -239,6 +240,16 @@ SUBROUTINE INITIALIZE(timeP) pSa = 0._RKIND END IF +! Variable wall properties - SCHWARZ July 2021---------------------- +! VARIABLE WALL + IF (useVarWall) THEN + IF (ALLOCATED(vWP0)) err = "Variable wall properties "// + 2 "already allocated. Correction needed" + ALLOCATE(vWP0(nvwp,tnNo)) + vWP0 = 0._RKIND + END IF +! ------------------------------------------------------------------ + ! Electrophysiology IF (cepEq) THEN ALLOCATE(Xion(nXion,tnNo)) @@ -663,6 +674,10 @@ SUBROUTINE FINALIZE IF (ALLOCATED(pSn)) DEALLOCATE(pSn) IF (ALLOCATED(pSa)) DEALLOCATE(pSa) +! Variable wall properties - SCHWARZ July 2021 --------------------- + IF (ALLOCATED(vWP0)) DEALLOCATE(vWP0) +! ------------------------------------------------------------------ + IF (ALLOCATED(Pinit)) DEALLOCATE(Pinit) IF (ALLOCATED(Vinit)) DEALLOCATE(Vinit) IF (ALLOCATED(Dinit)) DEALLOCATE(Dinit) diff --git a/Code/Source/svFSI/READFILES.f b/Code/Source/svFSI/READFILES.f index 0fa13395..22e3ba10 100755 --- a/Code/Source/svFSI/READFILES.f +++ b/Code/Source/svFSI/READFILES.f @@ -97,6 +97,9 @@ SUBROUTINE READFILES pstEq = .FALSE. sstEq = .FALSE. ibFlag = .FALSE. +! Variable wall properties - SCHWARZ------------------ + useVarWall = .FALSE. +! ---------------------------------------------------- i = IARGC() IF (i .NE. 0) THEN diff --git a/Code/Source/svFSI/REMESH.f b/Code/Source/svFSI/REMESH.f index 2911ab11..4bbd7bc1 100755 --- a/Code/Source/svFSI/REMESH.f +++ b/Code/Source/svFSI/REMESH.f @@ -394,6 +394,10 @@ END SUBROUTINE INTERP IF (ALLOCATED(pSn)) DEALLOCATE(pSn) IF (ALLOCATED(pSa)) DEALLOCATE(pSa) +! Variable wall properties - SCHWARZ July 2021 --------------------- + IF (ALLOCATED(vWP0)) DEALLOCATE(vWP0) +! ------------------------------------------------------------------ + t2 = CPUT() std = " Time taken for remeshing: "//STR(t2-t1)//" (s)" From c3da8c4b72e64c5e1239bcd69ff22bcf5c001146 Mon Sep 17 00:00:00 2001 From: Erica Date: Tue, 27 Jul 2021 12:40:01 -0700 Subject: [PATCH 03/10] Variable walls on multiple processors --- Code/Source/svFSI/DISTRIBUTE.f | 1 + Code/Source/svFSI/INITIALIZE.f | 12 +----------- Code/Source/svFSI/STRUCT.f | 25 +++++++++++++++++++------ 3 files changed, 21 insertions(+), 17 deletions(-) diff --git a/Code/Source/svFSI/DISTRIBUTE.f b/Code/Source/svFSI/DISTRIBUTE.f index a13c4e20..20b6a870 100755 --- a/Code/Source/svFSI/DISTRIBUTE.f +++ b/Code/Source/svFSI/DISTRIBUTE.f @@ -188,6 +188,7 @@ SUBROUTINE DISTRIBUTE CALL cm%bcast(pstEq) ! Variable wall properties - SCHWARZ July 2021------------------- CALL cm%bcast(useVarWall) + CALL cm%bcast(nvwp) ! --------------------------------------------------------------- CALL cm%bcast(sstEq) CALL cm%bcast(cepEq) diff --git a/Code/Source/svFSI/INITIALIZE.f b/Code/Source/svFSI/INITIALIZE.f index c35fb417..924f2d14 100755 --- a/Code/Source/svFSI/INITIALIZE.f +++ b/Code/Source/svFSI/INITIALIZE.f @@ -154,7 +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 (useVarWall) i = i + nvwp IF (cepEq) THEN i = i + nXion IF (cem%cpld) i = i + 1 @@ -240,16 +240,6 @@ SUBROUTINE INITIALIZE(timeP) pSa = 0._RKIND END IF -! Variable wall properties - SCHWARZ July 2021---------------------- -! VARIABLE WALL - IF (useVarWall) THEN - IF (ALLOCATED(vWP0)) err = "Variable wall properties "// - 2 "already allocated. Correction needed" - ALLOCATE(vWP0(nvwp,tnNo)) - vWP0 = 0._RKIND - END IF -! ------------------------------------------------------------------ - ! Electrophysiology IF (cepEq) THEN ALLOCATE(Xion(nXion,tnNo)) diff --git a/Code/Source/svFSI/STRUCT.f b/Code/Source/svFSI/STRUCT.f index 0f6d7ac5..a3d8e397 100755 --- a/Code/Source/svFSI/STRUCT.f +++ b/Code/Source/svFSI/STRUCT.f @@ -64,6 +64,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) ! Loop over all elements of mesh DO e=1, lM%nEl + ! Update domain and proceed if domain phys and eqn phys match cDmn = DOMAIN(lM, cEq, e) cPhys = eq(cEq)%dmn(cDmn)%phys @@ -87,9 +88,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) ! Variable wall - SCHWARZ July 2021--------------------------- ! Calculate local wall property - IF (useVarWall) THEN - lVWP(:,a) = vWP0(:,Ac) - END IF + IF (useVarWall) lVWP(:,a) = vWP0(:,Ac) ! ------------------------------------------------------------ IF (ALLOCATED(lM%fN)) THEN DO iFn=1, nFn @@ -190,6 +189,11 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, ya_g = 0._RKIND eVWP = 0._RKIND +! PRINT *, "Beginning of loop" +! PRINT *, eVWP(1) +! PRINT *, eVWP(2) +! PRINT *, eVWP(3) + DO a=1, eNoN ud(1) = ud(1) + N(a)*(rho*(al(i,a)-bfl(1,a)) + dmp*yl(i,a)) ud(2) = ud(2) + N(a)*(rho*(al(j,a)-bfl(2,a)) + dmp*yl(j,a)) @@ -207,9 +211,12 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, ! Variable wall - SCHWARZ July 2021--------------------------------- ! Calculate local wall property - IF (useVarWall) THEN - eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) - END IF + IF (useVarWall) eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) + +! PRINT *, "In loop" +! PRINT *, eVWP(1) +! PRINT *, eVWP(2) +! PRINT *, eVWP(3) ! ------------------------------------------------------------------ S0(1,1) = S0(1,1) + N(a)*pS0l(1,a) @@ -225,6 +232,12 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, S0(3,2) = S0(2,3) S0(1,3) = S0(3,1) +! PRINT *, "End of loop" +! PRINT *, eVWP(1) +! PRINT *, eVWP(2) +! PRINT *, eVWP(3) + + ! 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in ! Voigt notationa (Dm) CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm, eVWP) From 5c4ec175f7484a75714ac016d1684f85182fbf25 Mon Sep 17 00:00:00 2001 From: Erica Date: Wed, 4 Aug 2021 11:18:11 -0700 Subject: [PATCH 04/10] Adding LELAS varwall --- Code/Source/svFSI/FSI.f | 14 ++++++--- Code/Source/svFSI/LELAS.f | 60 +++++++++++++++++++++++++++++-------- Code/Source/svFSI/READMSH.f | 1 + Code/Source/svFSI/STRUCT.f | 17 +---------- 4 files changed, 59 insertions(+), 33 deletions(-) diff --git a/Code/Source/svFSI/FSI.f b/Code/Source/svFSI/FSI.f index 56bc3211..4261595f 100755 --- a/Code/Source/svFSI/FSI.f +++ b/Code/Source/svFSI/FSI.f @@ -50,7 +50,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(:), N(:), - 3 Nx(:,:), lR(:,:), lK(:,:,:), lKd(:,:,:) + 3 Nx(:,:), lR(:,:), lK(:,:,:), lKd(:,:,:), lVWP(:,:) eNoN = lM%eNoN nFn = lM%nFn @@ -59,7 +59,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), N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), - 3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN)) + 4 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN), lVWP(nvwp,eNoN)) ! Loop over all elements of mesh DO e=1, lM%nEl @@ -85,6 +85,12 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg) yl(:,a) = Yg(:,Ac) dl(:,a) = Dg(:,Ac) bfl(:,a) = Bf(:,Ac) + +! Variable wall - SCHWARZ July 2021--------------------------- +! 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) @@ -118,11 +124,11 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg) CASE (phys_lElas) CALL LELAS3D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, - 2 lR, lK) + 2 lR, lK, lVWP) CASE (phys_struct) CALL STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, - 2 fN, pS0l, pSl, ya_l, lR, lK) + 2 fN, pS0l, pSl, ya_l, lR, lK, lVWP) CASE (phys_ustruct) c CALL USTRUCT3D(eNoN, nFn, w, Jac, N, Nx, al, yl, dl, diff --git a/Code/Source/svFSI/LELAS.f b/Code/Source/svFSI/LELAS.f index c2a15013..bea1f7b8 100755 --- a/Code/Source/svFSI/LELAS.f +++ b/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) +! Variable wall - SCHWARZ July 2021--------------------------- +! 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,18 +138,19 @@ 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) @@ -161,6 +172,8 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, 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,6 +186,11 @@ 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) +! Variable wall - SCHWARZ July 2021--------------------------------- +! Calculate local wall property + IF (useVarWall) eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) +! ------------------------------------------------------------------ + 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) @@ -182,13 +200,29 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, END DO 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 + +! Variable wall - SCHWARZ July 2021--------------------------------- +! Calculate local wall property + IF (useVarWall) THEN + Cst(1,:) = eVWP(1:6) + Cst(2,:) = eVWP(7:12) + Cst(3,:) = eVWP(13:18) + Cst(4,:) = eVWP(19:24) + Cst(5,:) = eVWP(25:30) + Cst(6,:) = eVWP(31:36) + S(:) = MATMUL(Cst,ed) + ELSE + ! 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 + END IF + +! ------------------------------------------------------------------ + pSl = S ! Add prestress contribution diff --git a/Code/Source/svFSI/READMSH.f b/Code/Source/svFSI/READMSH.f index c0177cc2..b6fd1462 100755 --- a/Code/Source/svFSI/READMSH.f +++ b/Code/Source/svFSI/READMSH.f @@ -421,6 +421,7 @@ SUBROUTINE READMSH(list) END IF flag = .TRUE. useVarWall = .TRUE. + PRINT *,"Setting varwall flag to true" ALLOCATE(msh(iM)%x(nvwp,msh(iM)%gnNo)) msh(iM)%x = 0._RKIND CALL READVTUPDATA(msh(iM), cTmp, "varWallProps", nvwp, 1) diff --git a/Code/Source/svFSI/STRUCT.f b/Code/Source/svFSI/STRUCT.f index a3d8e397..19d9aba4 100755 --- a/Code/Source/svFSI/STRUCT.f +++ b/Code/Source/svFSI/STRUCT.f @@ -62,6 +62,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) 3 pSl(nsymd), ya_l(eNoN), N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), 4 lK(dof*dof,eNoN,eNoN), lVWP(nvwp,eNoN)) + ! Loop over all elements of mesh DO e=1, lM%nEl @@ -189,11 +190,6 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, ya_g = 0._RKIND eVWP = 0._RKIND -! PRINT *, "Beginning of loop" -! PRINT *, eVWP(1) -! PRINT *, eVWP(2) -! PRINT *, eVWP(3) - DO a=1, eNoN ud(1) = ud(1) + N(a)*(rho*(al(i,a)-bfl(1,a)) + dmp*yl(i,a)) ud(2) = ud(2) + N(a)*(rho*(al(j,a)-bfl(2,a)) + dmp*yl(j,a)) @@ -212,11 +208,6 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, ! Variable wall - SCHWARZ July 2021--------------------------------- ! Calculate local wall property IF (useVarWall) eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) - -! PRINT *, "In loop" -! PRINT *, eVWP(1) -! PRINT *, eVWP(2) -! PRINT *, eVWP(3) ! ------------------------------------------------------------------ S0(1,1) = S0(1,1) + N(a)*pS0l(1,a) @@ -232,12 +223,6 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, S0(3,2) = S0(2,3) S0(1,3) = S0(3,1) -! PRINT *, "End of loop" -! PRINT *, eVWP(1) -! PRINT *, eVWP(2) -! PRINT *, eVWP(3) - - ! 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in ! Voigt notationa (Dm) CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm, eVWP) From 9193621498198d557a441a50df48a89ade7634e0 Mon Sep 17 00:00:00 2001 From: Erica Date: Wed, 4 Aug 2021 15:23:07 -0700 Subject: [PATCH 05/10] Bug fixing for varwall in MESH --- Code/Source/svFSI/EQASSEM.f | 2 ++ Code/Source/svFSI/MATMODELS.f | 4 ---- Code/Source/svFSI/MESH.f | 11 ++++++++--- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/Code/Source/svFSI/EQASSEM.f b/Code/Source/svFSI/EQASSEM.f index 49ac9408..8d25f6d6 100755 --- a/Code/Source/svFSI/EQASSEM.f +++ b/Code/Source/svFSI/EQASSEM.f @@ -42,6 +42,8 @@ SUBROUTINE GLOBALEQASSEM(lM, Ag, Yg, Dg) REAL(KIND=RKIND), INTENT(IN) :: Ag(tDof,tnNo), Yg(tDof,tnNo), 2 Dg(tDof,tnNo) + PRINT *, eq(cEq)%phys + SELECT CASE (eq(cEq)%phys) CASE (phys_fluid) CALL CONSTRUCT_FLUID(lM, Ag, Yg) diff --git a/Code/Source/svFSI/MATMODELS.f b/Code/Source/svFSI/MATMODELS.f index 0eaf7497..712b7e2b 100755 --- a/Code/Source/svFSI/MATMODELS.f +++ b/Code/Source/svFSI/MATMODELS.f @@ -135,12 +135,8 @@ SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm, eVWP) IF (useVarWall) THEN ! Converting elastic modulus and poisson ratio to g1 g1 = eVWP(1)* 0.5_RKIND/(1._RKIND+eVWP(2)) -! PRINT *, "USING VARWALL" -! PRINT *, g1 ELSE g1 = 2._RKIND * stM%C10 -! PRINT *, "NOT USING VARWALL" -! PRINT *, g1 END IF Sb = g1*IDm diff --git a/Code/Source/svFSI/MESH.f b/Code/Source/svFSI/MESH.f index e77ddd71..3a5cdbca 100755 --- a/Code/Source/svFSI/MESH.f +++ b/Code/Source/svFSI/MESH.f @@ -49,7 +49,7 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:) REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), dl(:,:), 2 dol(:,:), bfl(:,:), pS0l(:,:), pSl(:), N(:), Nx(:,:), lR(:,:), - 3 lK(:,:,:) + 3 lK(:,:,:), lVWP(:,:) eNoN = lM%eNoN IF (.NOT.mvMsh) err = "Mesh equation is solved for moving mesh/"// @@ -58,7 +58,8 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) ! MESH: dof = nsd ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), dl(tDof,eNoN), 2 dol(nsd,eNoN), bfl(nsd,eNoN), pS0l(nsymd,eNoN), pSl(nsymd), - 3 N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN)) + 3 N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN), + 4 lVWP(nvwp,eNoN)) ! Start and end DOF is = nsd + 2 @@ -84,6 +85,10 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) al(:,a) = Ag(:,Ac) dl(:,a) = Dg(:,Ac) dol(:,a) = Do(is:ie,Ac) +! Variable wall - SCHWARZ July 2021--------------------------- +! Calculate local wall property + IF (useVarWall) lVWP(:,a) = vWP0(:,Ac) +! ------------------------------------------------------------ END DO ! For MESH, the reference configuration is the one at the @@ -105,7 +110,7 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) pS0l = 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, From b76fc12fc3217cc31767013109f88704748af6e7 Mon Sep 17 00:00:00 2001 From: Erica Date: Wed, 4 Aug 2021 16:57:03 -0700 Subject: [PATCH 06/10] Fixed MESH exception to using varwall properties --- Code/Source/svFSI/EQASSEM.f | 2 -- Code/Source/svFSI/LELAS.f | 7 +++++-- Code/Source/svFSI/MESH.f | 11 +++-------- 3 files changed, 8 insertions(+), 12 deletions(-) diff --git a/Code/Source/svFSI/EQASSEM.f b/Code/Source/svFSI/EQASSEM.f index 8d25f6d6..49ac9408 100755 --- a/Code/Source/svFSI/EQASSEM.f +++ b/Code/Source/svFSI/EQASSEM.f @@ -42,8 +42,6 @@ SUBROUTINE GLOBALEQASSEM(lM, Ag, Yg, Dg) REAL(KIND=RKIND), INTENT(IN) :: Ag(tDof,tnNo), Yg(tDof,tnNo), 2 Dg(tDof,tnNo) - PRINT *, eq(cEq)%phys - SELECT CASE (eq(cEq)%phys) CASE (phys_fluid) CALL CONSTRUCT_FLUID(lM, Ag, Yg) diff --git a/Code/Source/svFSI/LELAS.f b/Code/Source/svFSI/LELAS.f index bea1f7b8..cc6023fb 100755 --- a/Code/Source/svFSI/LELAS.f +++ b/Code/Source/svFSI/LELAS.f @@ -188,7 +188,10 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, ! Variable wall - SCHWARZ July 2021--------------------------------- ! Calculate local wall property - IF (useVarWall) eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) +! Don't use on the mesh part though lol + 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) @@ -203,7 +206,7 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, ! Variable wall - SCHWARZ July 2021--------------------------------- ! Calculate local wall property - IF (useVarWall) THEN + IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN Cst(1,:) = eVWP(1:6) Cst(2,:) = eVWP(7:12) Cst(3,:) = eVWP(13:18) diff --git a/Code/Source/svFSI/MESH.f b/Code/Source/svFSI/MESH.f index 3a5cdbca..e77ddd71 100755 --- a/Code/Source/svFSI/MESH.f +++ b/Code/Source/svFSI/MESH.f @@ -49,7 +49,7 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:) REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), dl(:,:), 2 dol(:,:), bfl(:,:), pS0l(:,:), pSl(:), N(:), Nx(:,:), lR(:,:), - 3 lK(:,:,:), lVWP(:,:) + 3 lK(:,:,:) eNoN = lM%eNoN IF (.NOT.mvMsh) err = "Mesh equation is solved for moving mesh/"// @@ -58,8 +58,7 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) ! MESH: dof = nsd ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), dl(tDof,eNoN), 2 dol(nsd,eNoN), bfl(nsd,eNoN), pS0l(nsymd,eNoN), pSl(nsymd), - 3 N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN), - 4 lVWP(nvwp,eNoN)) + 3 N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN), lK(dof*dof,eNoN,eNoN)) ! Start and end DOF is = nsd + 2 @@ -85,10 +84,6 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) al(:,a) = Ag(:,Ac) dl(:,a) = Dg(:,Ac) dol(:,a) = Do(is:ie,Ac) -! Variable wall - SCHWARZ July 2021--------------------------- -! Calculate local wall property - IF (useVarWall) lVWP(:,a) = vWP0(:,Ac) -! ------------------------------------------------------------ END DO ! For MESH, the reference configuration is the one at the @@ -110,7 +105,7 @@ SUBROUTINE CONSTRUCT_MESH(lM, Ag, Dg) pS0l = 0._RKIND IF (nsd .EQ. 3) THEN CALL LELAS3D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, lR, - 2 lK,lVWP) + 2 lK) ELSE IF (nsd .EQ. 2) THEN CALL LELAS2D(eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, lR, From 82f86b7589cc0c220b4b62e1c77ddfb9682c3116 Mon Sep 17 00:00:00 2001 From: Erica Date: Sun, 8 Aug 2021 18:58:57 -0700 Subject: [PATCH 07/10] Added several input checks for variable wall ALE --- Code/Source/svFSI/POST.f | 50 ++++++++++++++++++++++++++--------- Code/Source/svFSI/READFILES.f | 24 +++++++++++++++++ 2 files changed, 61 insertions(+), 13 deletions(-) diff --git a/Code/Source/svFSI/POST.f b/Code/Source/svFSI/POST.f index e89d45e9..78f08941 100755 --- a/Code/Source/svFSI/POST.f +++ b/Code/Source/svFSI/POST.f @@ -521,7 +521,8 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) REAL(KIND=RKIND) w, Jac, detF, Je, ya, Ja, elM, nu, lambda, mu, 2 p, trS, vmises, xi(nsd), xi0(nsd), xp(nsd), ed(nsymd), 3 Im(nsd,nsd), F(nsd,nsd), C(nsd,nsd), Eg(nsd,nsd), P1(nsd,nsd), - 4 S(nsd,nsd), sigma(nsd,nsd), Dm(nsymd,nsymd), eVWP(nvwp) + 4 S(nsd,nsd), sigma(nsd,nsd), Dm(nsymd,nsymd), eVWP(nvwp), + 5 sigma_temp(6), Cst(6,6) TYPE(fsType) :: fs INTEGER, ALLOCATABLE :: eNds(:) @@ -563,6 +564,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) sE = 0._RKIND insd = nsd ya = 0._RKIND + IF (lM%lFib) insd = 1 DO e=1, lM%nEl @@ -709,20 +711,42 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) CASE (outGrp_stress, outGrp_cauchy, outGrp_mises) sigma = 0._RKIND + sigma_temp = 0._RKIND IF (cPhys .EQ. phys_lElas) THEN IF (nsd .EQ. 3) THEN - detF = lambda*(ed(1) + ed(2) + ed(3)) - sigma(1,1) = detF + 2._RKIND*mu*ed(1) - sigma(2,2) = detF + 2._RKIND*mu*ed(2) - sigma(3,3) = detF + 2._RKIND*mu*ed(3) - - sigma(1,2) = mu*ed(4) - sigma(2,3) = mu*ed(5) - sigma(3,1) = mu*ed(6) - - sigma(2,1) = sigma(1,2) - sigma(3,2) = sigma(2,3) - sigma(1,3) = sigma(3,1) + IF (useVarWall) THEN + Cst(1,:) = eVWP(1:6) + Cst(2,:) = eVWP(7:12) + Cst(3,:) = eVWP(13:18) + Cst(4,:) = eVWP(19:24) + Cst(5,:) = eVWP(25:30) + Cst(6,:) = eVWP(31:36) + sigma_temp(:) = MATMUL(Cst,ed) + sigma(1,1) = sigma_temp(1) + sigma(2,2) = sigma_temp(2) + sigma(3,3) = sigma_temp(3) + + sigma(1,2) = sigma_temp(4) + sigma(2,3) = sigma_temp(5) + sigma(3,1) = sigma_temp(6) + + sigma(2,1) = sigma(1,2) + sigma(3,2) = sigma(2,3) + sigma(1,3) = sigma(3,1) + ELSE + detF = lambda*(ed(1) + ed(2) + ed(3)) + sigma(1,1) = detF + 2._RKIND*mu*ed(1) + sigma(2,2) = detF + 2._RKIND*mu*ed(2) + sigma(3,3) = detF + 2._RKIND*mu*ed(3) + + sigma(1,2) = mu*ed(4) + sigma(2,3) = mu*ed(5) + sigma(3,1) = mu*ed(6) + + sigma(2,1) = sigma(1,2) + sigma(3,2) = sigma(2,3) + sigma(1,3) = sigma(3,1) + END IF ELSE detF = lambda*(ed(1) + ed(2)) sigma(1,1) = detF + 2._RKIND*mu*ed(1) diff --git a/Code/Source/svFSI/READFILES.f b/Code/Source/svFSI/READFILES.f index 6617c380..2bc33593 100755 --- a/Code/Source/svFSI/READFILES.f +++ b/Code/Source/svFSI/READFILES.f @@ -417,6 +417,7 @@ SUBROUTINE READEQ(lEq, list, eqName) THflag = .FALSE. propL = prop_NA + SELECT CASE (eqName) ! FLUID Navier-Stokes solver ------------------------------------ CASE ('fluid') @@ -487,6 +488,11 @@ SUBROUTINE READEQ(lEq, list, eqName) IF (nsd .EQ. 3) propL(6,1) = f_z CALL READDOMAIN(lEq, propL, list) + IF (useVarWall .AND. (nvwp .LT. 36._RKIND)) THEN + err = "Number of variable wall properties for linear "// + 2 "elastic material must be at least 36." + END IF + lPtr => list%get(pstEq, "Prestress") IF (pstEq) THEN @@ -563,6 +569,8 @@ SUBROUTINE READEQ(lEq, list, eqName) IF (nsd .EQ. 3) propL(8,1) = f_z CALL READDOMAIN(lEq, propL, list) + IF (useVarWall) err = "Varible wall for USTRUCT is not "// + 2 "implemented yet" lPtr => list%get(pstEq, "Prestress") IF (pstEq) err = "Prestress for USTRUCT is not implemented yet" @@ -601,6 +609,9 @@ SUBROUTINE READEQ(lEq, list, eqName) propL(8,1) = f_z CALL READDOMAIN(lEq, propL, list) + IF (useVarWall) err = "Varible wall for SHELLS is not "// + 2 "implemented yet" + lPtr => list%get(pstEq, "Prestress") IF (pstEq) err = "Prestress for SHELLS is not implemented yet" @@ -613,6 +624,7 @@ SUBROUTINE READEQ(lEq, list, eqName) ! COUPLED MOMENTUM FLUID STRUCTURE INTERACTION equation solver--- CASE ('CMM') + IF (nsd .NE. 3) err = "CMM eq. is not implemented for 2D "// 2 "domains" lEq%phys = phys_CMM @@ -704,6 +716,9 @@ SUBROUTINE READEQ(lEq, list, eqName) CALL READDOMAIN(lEq, propL, list) IF (cmmInit) lEq%dmn(:)%prop(solid_density) = 0._RKIND + IF (useVarWall) err = "Varible wall for CMM should not be "// + 2 "defined in mesh data block" + CALL READLS(lSolver_GMRES, lEq, list) ! FLUID STRUCTURE INTERACTION equation solver--------------------- @@ -2490,6 +2505,15 @@ SUBROUTINE READMATMODEL(lDmn, lPD) err = "Undefined constitutive model used" END SELECT + IF (useVarWall .AND. (lDmn%stM%isoType .NE. stIso_nHook)) THEN + err = "Variable wall properties currently only implemented "// + 2 "for isotropic Neo-Hookean material in STRUCT." + END IF + IF (useVarWall .AND. (nvwp .LT. 2._RKIND)) THEN + err = "Number of variable wall properties for isotropic "// + 2 "Neo-Hookean must be at least 2." + END IF + ! Fiber reinforcement stress lFib => lPD%get(ctmp, "Fiber reinforcement stress") IF (ASSOCIATED(lFib)) THEN From 21206ab423bd65afd825685af0f497a2ab91995b Mon Sep 17 00:00:00 2001 From: Erica Date: Sun, 8 Aug 2021 21:06:07 -0700 Subject: [PATCH 08/10] Added cauchy output to LELAS material- for some reason sigma was assigned but not allowed to be outputed (may have been intended to be assigned to Stress as S variable, but that is never used with LELAS) --- Code/Source/svFSI/POST.f | 2 -- Code/Source/svFSI/READFILES.f | 3 ++- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/Code/Source/svFSI/POST.f b/Code/Source/svFSI/POST.f index 78f08941..28d5502f 100755 --- a/Code/Source/svFSI/POST.f +++ b/Code/Source/svFSI/POST.f @@ -710,8 +710,6 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) END IF CASE (outGrp_stress, outGrp_cauchy, outGrp_mises) - sigma = 0._RKIND - sigma_temp = 0._RKIND IF (cPhys .EQ. phys_lElas) THEN IF (nsd .EQ. 3) THEN IF (useVarWall) THEN diff --git a/Code/Source/svFSI/READFILES.f b/Code/Source/svFSI/READFILES.f index 2bc33593..bcd22e2a 100755 --- a/Code/Source/svFSI/READFILES.f +++ b/Code/Source/svFSI/READFILES.f @@ -501,7 +501,7 @@ SUBROUTINE READEQ(lEq, list, eqName) outPuts(2) = out_stress outPuts(3) = out_strain ELSE - nDOP = (/8,2,0,0/) + nDOP = (/9,2,0,0/) outPuts(1) = out_displacement outPuts(2) = out_mises outPuts(3) = out_stress @@ -510,6 +510,7 @@ SUBROUTINE READEQ(lEq, list, eqName) outPuts(6) = out_acceleration outPuts(7) = out_integ outPuts(8) = out_jacobian + outPuts(9) = out_cauchy END IF CALL READLS(lSolver_CG, lEq, list) From c6adf3e7800ec4d42045ea5997b2cd04134df398 Mon Sep 17 00:00:00 2001 From: Erica Date: Wed, 15 Sep 2021 11:21:08 -0700 Subject: [PATCH 09/10] Finalized isotropic input and post-processing --- Code/Source/svFSI/LELAS.f | 62 ++++++++++++++++------------------- Code/Source/svFSI/MOD.f | 2 +- Code/Source/svFSI/POST.f | 54 +++++++++++------------------- Code/Source/svFSI/READFILES.f | 34 ++++++++++++------- 4 files changed, 73 insertions(+), 79 deletions(-) diff --git a/Code/Source/svFSI/LELAS.f b/Code/Source/svFSI/LELAS.f index cc6023fb..013d4572 100755 --- a/Code/Source/svFSI/LELAS.f +++ b/Code/Source/svFSI/LELAS.f @@ -153,8 +153,6 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, 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) @@ -162,13 +160,6 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, 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 @@ -186,9 +177,9 @@ 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) -! Variable wall - SCHWARZ July 2021--------------------------------- +! ------------------------------------------------------------------ ! Calculate local wall property -! Don't use on the mesh part though lol +! Don't use if calculating mesh IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) END IF @@ -201,30 +192,33 @@ PURE SUBROUTINE LELAS3D (eNoN, w, N, Nx, al, dl, bfl, pS0l, pSl, S0(5) = S0(5) + N(a)*pS0l(5,a) S0(6) = S0(6) + N(a)*pS0l(6,a) END DO - divD = lambda*(ed(1) + ed(2) + ed(3)) -! Variable wall - SCHWARZ July 2021--------------------------------- -! Calculate local wall property IF (useVarWall .AND. (phys_mesh .NE. eq(cEq)%phys)) THEN - Cst(1,:) = eVWP(1:6) - Cst(2,:) = eVWP(7:12) - Cst(3,:) = eVWP(13:18) - Cst(4,:) = eVWP(19:24) - Cst(5,:) = eVWP(25:30) - Cst(6,:) = eVWP(31:36) - S(:) = MATMUL(Cst,ed) + elM = eVWP(1) + nu = evWP(2) ELSE - ! 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 + 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 @@ -290,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) @@ -333,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)) @@ -343,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) @@ -354,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 diff --git a/Code/Source/svFSI/MOD.f b/Code/Source/svFSI/MOD.f index ebdc78be..d85bce66 100755 --- a/Code/Source/svFSI/MOD.f +++ b/Code/Source/svFSI/MOD.f @@ -941,7 +941,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(:,:) diff --git a/Code/Source/svFSI/POST.f b/Code/Source/svFSI/POST.f index 28d5502f..1b396baa 100755 --- a/Code/Source/svFSI/POST.f +++ b/Code/Source/svFSI/POST.f @@ -559,6 +559,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) 2 dl(tDof,fs%eNoN), yl(tDof,fs%eNoN), fN(nsd,nFn), resl(m), 3 Nx(nsd,fs%eNoN), N(fs%eNoN), lVWP(nvwp,fs%eNoN)) + S = 0._RKIND sA = 0._RKIND sF = 0._RKIND sE = 0._RKIND @@ -574,7 +575,7 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) 2 cPhys .NE. phys_ustruct .AND. 3 cPhys .NE. phys_lElas) CYCLE - IF (cPhys .EQ. phys_lElas) THEN + IF (cPhys .EQ. phys_lElas .AND. .NOT. useVarWall) THEN elM = eq(iEq)%dmn(cDmn)%prop(elasticity_modulus) nu = eq(iEq)%dmn(cDmn)%prop(poisson_ratio) lambda = elM*nu/(1._RKIND + nu)/(1._RKIND - 2._RKIND*nu) @@ -642,6 +643,12 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) F(2,2) = F(2,2) + Nx(2,a)*dl(j,a) END IF END DO + IF (cPhys .EQ. phys_lElas .AND. useVarWall) THEN + elM = eVWP(1) + nu = eVWP(2) + lambda = elM*nu/(1._RKIND + nu)/(1._RKIND - 2._RKIND*nu) + mu = 0.5_RKIND*elM/(1._RKIND + nu) + END IF detF = MAT_DET(F, nsd) ed = 0._RKIND @@ -712,39 +719,18 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) CASE (outGrp_stress, outGrp_cauchy, outGrp_mises) IF (cPhys .EQ. phys_lElas) THEN IF (nsd .EQ. 3) THEN - IF (useVarWall) THEN - Cst(1,:) = eVWP(1:6) - Cst(2,:) = eVWP(7:12) - Cst(3,:) = eVWP(13:18) - Cst(4,:) = eVWP(19:24) - Cst(5,:) = eVWP(25:30) - Cst(6,:) = eVWP(31:36) - sigma_temp(:) = MATMUL(Cst,ed) - sigma(1,1) = sigma_temp(1) - sigma(2,2) = sigma_temp(2) - sigma(3,3) = sigma_temp(3) - - sigma(1,2) = sigma_temp(4) - sigma(2,3) = sigma_temp(5) - sigma(3,1) = sigma_temp(6) - - sigma(2,1) = sigma(1,2) - sigma(3,2) = sigma(2,3) - sigma(1,3) = sigma(3,1) - ELSE - detF = lambda*(ed(1) + ed(2) + ed(3)) - sigma(1,1) = detF + 2._RKIND*mu*ed(1) - sigma(2,2) = detF + 2._RKIND*mu*ed(2) - sigma(3,3) = detF + 2._RKIND*mu*ed(3) - - sigma(1,2) = mu*ed(4) - sigma(2,3) = mu*ed(5) - sigma(3,1) = mu*ed(6) - - sigma(2,1) = sigma(1,2) - sigma(3,2) = sigma(2,3) - sigma(1,3) = sigma(3,1) - END IF + detF = lambda*(ed(1) + ed(2) + ed(3)) + sigma(1,1) = detF + 2._RKIND*mu*ed(1) + sigma(2,2) = detF + 2._RKIND*mu*ed(2) + sigma(3,3) = detF + 2._RKIND*mu*ed(3) + + sigma(1,2) = mu*ed(4) + sigma(2,3) = mu*ed(5) + sigma(3,1) = mu*ed(6) + + sigma(2,1) = sigma(1,2) + sigma(3,2) = sigma(2,3) + sigma(1,3) = sigma(3,1) ELSE detF = lambda*(ed(1) + ed(2)) sigma(1,1) = detF + 2._RKIND*mu*ed(1) diff --git a/Code/Source/svFSI/READFILES.f b/Code/Source/svFSI/READFILES.f index bcd22e2a..58c2a9c7 100755 --- a/Code/Source/svFSI/READFILES.f +++ b/Code/Source/svFSI/READFILES.f @@ -226,6 +226,10 @@ SUBROUTINE READFILES !-------------------------------------------------------------------- ! Reading equations + + IF ((nsd .NE. 3) .AND. useVarWall) err = "Variable wall "// + 2 "properties can only be used in 3-dimensions" + nEq = list%srch("Add equation",ll=1) std = " Number of equations: "//nEq ALLOCATE(eq(nEq)) @@ -431,7 +435,7 @@ SUBROUTINE READEQ(lEq, list, eqName) IF (nsd .EQ. 3) propL(5,1) = f_z CALL READDOMAIN(lEq, propL, list) - nDOP = (/11,2,3,0/) + nDOP = (/12,2,3,0/) outPuts(1) = out_velocity outPuts(2) = out_pressure outPuts(3) = out_WSS @@ -444,6 +448,7 @@ SUBROUTINE READEQ(lEq, list, eqName) outPuts(9) = out_viscosity outPuts(10) = out_divergence outPuts(11) = out_acceleration + outPuts(12) = out_displacement CALL READLS(lSolver_NS, lEq, list) @@ -480,19 +485,26 @@ SUBROUTINE READEQ(lEq, list, eqName) CASE ('lElas') lEq%phys = phys_lElas - propL(1,1) = solid_density - propL(2,1) = elasticity_modulus - propL(3,1) = poisson_ratio - propL(4,1) = f_x - propL(5,1) = f_y - IF (nsd .EQ. 3) propL(6,1) = f_z - CALL READDOMAIN(lEq, propL, list) - - IF (useVarWall .AND. (nvwp .LT. 36._RKIND)) THEN + IF (useVarWall .AND. (nvwp .LT. 2._RKIND)) THEN err = "Number of variable wall properties for linear "// - 2 "elastic material must be at least 36." + 2 "elastic material must be at least 2." END IF + IF (useVarWall) THEN + propL(1,1) = solid_density + propL(2,1) = f_x + propL(3,1) = f_y + IF (nsd .EQ. 3) propL(4,1) = f_z + ELSE + propL(1,1) = solid_density + propL(2,1) = elasticity_modulus + propL(3,1) = poisson_ratio + propL(4,1) = f_x + propL(5,1) = f_y + IF (nsd .EQ. 3) propL(6,1) = f_z + END IF + CALL READDOMAIN(lEq, propL, list) + lPtr => list%get(pstEq, "Prestress") IF (pstEq) THEN From 68469a02c861f72ea2a44c158a352cf33262499c Mon Sep 17 00:00:00 2001 From: Erica Date: Fri, 17 Sep 2021 18:30:08 -0700 Subject: [PATCH 10/10] Final cleanup on varwall inputs --- Code/Source/svFSI/DISTRIBUTE.f | 4 ++-- Code/Source/svFSI/FSI.f | 2 +- Code/Source/svFSI/INITIALIZE.f | 2 +- Code/Source/svFSI/LELAS.f | 2 +- Code/Source/svFSI/MATMODELS.f | 2 -- Code/Source/svFSI/MOD.f | 7 +------ Code/Source/svFSI/POST.f | 6 ++---- Code/Source/svFSI/READFILES.f | 2 -- Code/Source/svFSI/READMSH.f | 4 ++-- Code/Source/svFSI/REMESH.f | 2 +- Code/Source/svFSI/STRUCT.f | 4 ++-- 11 files changed, 13 insertions(+), 24 deletions(-) diff --git a/Code/Source/svFSI/DISTRIBUTE.f b/Code/Source/svFSI/DISTRIBUTE.f index 20b6a870..3e3a6de7 100755 --- a/Code/Source/svFSI/DISTRIBUTE.f +++ b/Code/Source/svFSI/DISTRIBUTE.f @@ -186,7 +186,7 @@ SUBROUTINE DISTRIBUTE CALL cm%bcast(cmmVarWall) CALL cm%bcast(shlEq) CALL cm%bcast(pstEq) -! Variable wall properties - SCHWARZ July 2021------------------- +! Varwall properties--------------------------------------------- CALL cm%bcast(useVarWall) CALL cm%bcast(nvwp) ! --------------------------------------------------------------- @@ -260,7 +260,7 @@ SUBROUTINE DISTRIBUTE DEALLOCATE(tmpX) END IF -! Variable wall properties - SCHWARZ July 2021 --------------------- +! Varwall properties------------------------------------------------ ! Distribute variable wall properties (vWP0) to processors flag = ALLOCATED(vWP0) CALL cm%bcast(flag) diff --git a/Code/Source/svFSI/FSI.f b/Code/Source/svFSI/FSI.f index 8bc9ca49..56a88f8b 100755 --- a/Code/Source/svFSI/FSI.f +++ b/Code/Source/svFSI/FSI.f @@ -99,7 +99,7 @@ SUBROUTINE CONSTRUCT_FSI(lM, Ag, Yg, Dg) dl(:,a) = Dg(:,Ac) bfl(:,a) = Bf(:,Ac) -! Variable wall - SCHWARZ July 2021--------------------------- +! Varwall properties ----------------------------------------- ! Calculate local wall property IF (useVarWall) lVWP(:,a) = vWP0(:,Ac) ! ------------------------------------------------------------ diff --git a/Code/Source/svFSI/INITIALIZE.f b/Code/Source/svFSI/INITIALIZE.f index 924f2d14..9379cfe3 100755 --- a/Code/Source/svFSI/INITIALIZE.f +++ b/Code/Source/svFSI/INITIALIZE.f @@ -664,7 +664,7 @@ SUBROUTINE FINALIZE IF (ALLOCATED(pSn)) DEALLOCATE(pSn) IF (ALLOCATED(pSa)) DEALLOCATE(pSa) -! Variable wall properties - SCHWARZ July 2021 --------------------- +! Varwall properties ----------------------------------------------- IF (ALLOCATED(vWP0)) DEALLOCATE(vWP0) ! ------------------------------------------------------------------ diff --git a/Code/Source/svFSI/LELAS.f b/Code/Source/svFSI/LELAS.f index 013d4572..940003db 100755 --- a/Code/Source/svFSI/LELAS.f +++ b/Code/Source/svFSI/LELAS.f @@ -79,7 +79,7 @@ SUBROUTINE CONSTRUCT_LELAS(lM, Ag, Dg) dl(:,a) = Dg(:,Ac) bfl(:,a) = Bf(:,Ac) IF (ALLOCATED(pS0)) pS0l(:,a) = pS0(:,Ac) -! Variable wall - SCHWARZ July 2021--------------------------- +! Varwall properties------------------------------------------ ! Calculate local wall property IF (useVarWall) THEN lVWP(:,a) = vWP0(:,Ac) diff --git a/Code/Source/svFSI/MATMODELS.f b/Code/Source/svFSI/MATMODELS.f index 712b7e2b..0842716a 100755 --- a/Code/Source/svFSI/MATMODELS.f +++ b/Code/Source/svFSI/MATMODELS.f @@ -46,9 +46,7 @@ SUBROUTINE GETPK2CC(lDmn, F, nfd, fl, ya, S, Dm, eVWP) 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) -! VARIABLE WALL PROPERTIES - SCHWARZ JULY 2021 --------------------- REAL(KIND=RKIND), INTENT(IN), OPTIONAL :: eVWP(nvwp) -! ------------------------------------------------------------------ TYPE(stModelType) :: stM REAL(KIND=RKIND) :: nd, Kp, J, J2d, J4d, trE, p, pl, Inv1, Inv2, diff --git a/Code/Source/svFSI/MOD.f b/Code/Source/svFSI/MOD.f index d85bce66..ccb4c0a6 100755 --- a/Code/Source/svFSI/MOD.f +++ b/Code/Source/svFSI/MOD.f @@ -837,10 +837,8 @@ MODULE COMMOD LOGICAL cmmInit ! Whether variable wall properties are used for CMM LOGICAL cmmVarWall -! SCHWARZ July 2021---------------------------------------------- ! Whether variable wall properties should be read and used LOGICAL useVarWall -! --------------------------------------------------------------- ! Whether shell equation is being solved LOGICAL shlEq ! Whether PRESTRESS is being solved @@ -897,10 +895,8 @@ MODULE COMMOD INTEGER(KIND=IKIND) rsTS ! Number of stress values to be stored INTEGER(KIND=IKIND) nsymd -! SCHWARZ July 2021---------------------------------------------- ! Number of variable wall properties to read in from mesh INTEGER(KIND=IKIND) nvwp -! --------------------------------------------------------------- ! REAL VARIABLES ! Time step size @@ -972,9 +968,8 @@ MODULE COMMOD REAL(KIND=RKIND), ALLOCATABLE :: pSn(:,:) REAL(KIND=RKIND), ALLOCATABLE :: pSa(:) -! Variables for variable wall properties - SCHWARZ July 2021---- +! Variables for variable wall properties REAL(KIND=RKIND), ALLOCATABLE :: vWP0(:,:) -! -------------------------------------------------------------- ! Temporary storage for initializing state variables REAL(KIND=RKIND), ALLOCATABLE :: Pinit(:) diff --git a/Code/Source/svFSI/POST.f b/Code/Source/svFSI/POST.f index 1b396baa..81477042 100755 --- a/Code/Source/svFSI/POST.f +++ b/Code/Source/svFSI/POST.f @@ -598,12 +598,11 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) xl(:,a) = x(:,Ac) dl(:,a) = lD(:,Ac) yl(:,a) = lY(:,Ac) -! Variable wall - SCHWARZ July 2021--------------------------- +! Varwall properties ----------------------------------------- ! Calculate local wall property IF (useVarWall) THEN lVWP(:,a) = vWP0(:,Ac) END IF -! ------------------------------------------------------------ END DO Je = 0._RKIND @@ -629,12 +628,11 @@ SUBROUTINE TPOST(lM, m, res, resE, lD, lY, iEq, outGrp) F(3,1) = F(3,1) + Nx(1,a)*dl(k,a) F(3,2) = F(3,2) + Nx(2,a)*dl(k,a) F(3,3) = F(3,3) + Nx(3,a)*dl(k,a) -! Variable wall - SCHWARZ July 2021--------------------- +! Variable wall ---------------------------------------- ! Calculate local wall property IF (useVarWall) THEN eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) END IF -! ------------------------------------------------------ ELSE F(1,1) = F(1,1) + Nx(1,a)*dl(i,a) diff --git a/Code/Source/svFSI/READFILES.f b/Code/Source/svFSI/READFILES.f index 58c2a9c7..77e1ddd5 100755 --- a/Code/Source/svFSI/READFILES.f +++ b/Code/Source/svFSI/READFILES.f @@ -97,9 +97,7 @@ SUBROUTINE READFILES pstEq = .FALSE. sstEq = .FALSE. ibFlag = .FALSE. -! Variable wall properties - SCHWARZ------------------ useVarWall = .FALSE. -! ---------------------------------------------------- i = IARGC() IF (i .NE. 0) THEN diff --git a/Code/Source/svFSI/READMSH.f b/Code/Source/svFSI/READMSH.f index 27c923d9..6220479a 100755 --- a/Code/Source/svFSI/READMSH.f +++ b/Code/Source/svFSI/READMSH.f @@ -408,7 +408,7 @@ SUBROUTINE READMSH(list) END DO END IF -! Read variable wall properties - SCHWARZ July 2021 ---------------- +! Read varwall properties ------------------------------------------ flag = .FALSE. DO iM=1, nMsh lPM => list%get(msh(iM)%name,"Add mesh",iM) @@ -421,7 +421,7 @@ SUBROUTINE READMSH(list) END IF flag = .TRUE. useVarWall = .TRUE. - PRINT *,"Setting varwall flag to true" + std = "Setting varwall flag to true" ALLOCATE(msh(iM)%x(nvwp,msh(iM)%gnNo)) msh(iM)%x = 0._RKIND CALL READVTUPDATA(msh(iM), cTmp, "varWallProps", nvwp, 1) diff --git a/Code/Source/svFSI/REMESH.f b/Code/Source/svFSI/REMESH.f index 4bbd7bc1..f693f0d7 100755 --- a/Code/Source/svFSI/REMESH.f +++ b/Code/Source/svFSI/REMESH.f @@ -394,7 +394,7 @@ END SUBROUTINE INTERP IF (ALLOCATED(pSn)) DEALLOCATE(pSn) IF (ALLOCATED(pSa)) DEALLOCATE(pSa) -! Variable wall properties - SCHWARZ July 2021 --------------------- +! Varwall properties ----------------------------------------------- IF (ALLOCATED(vWP0)) DEALLOCATE(vWP0) ! ------------------------------------------------------------------ diff --git a/Code/Source/svFSI/STRUCT.f b/Code/Source/svFSI/STRUCT.f index 19d9aba4..1afe9a67 100755 --- a/Code/Source/svFSI/STRUCT.f +++ b/Code/Source/svFSI/STRUCT.f @@ -87,7 +87,7 @@ SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg) dl(:,a) = Dg(:,Ac) bfl(:,a) = Bf(:,Ac) -! Variable wall - SCHWARZ July 2021--------------------------- +! Varwall properties----------------------------------------------------- ! Calculate local wall property IF (useVarWall) lVWP(:,a) = vWP0(:,Ac) ! ------------------------------------------------------------ @@ -205,7 +205,7 @@ SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, F(3,2) = F(3,2) + Nx(2,a)*dl(k,a) F(3,3) = F(3,3) + Nx(3,a)*dl(k,a) -! Variable wall - SCHWARZ July 2021--------------------------------- +! Varwall properits------------------------------------------------- ! Calculate local wall property IF (useVarWall) eVWP(:) = eVWP(:) + N(a)*lVWP(:,a) ! ------------------------------------------------------------------