from vb2py.vbfunctions import *
from vb2py.vbdebug import *
"""
CCYL1.bas : Nonlinear Concentric Cylinders : 1/2 Hybrid Scheme
Determine RLR Given Target Radial Cauchy Stress
Print Soution At Time Tau; Capture Shock Front Info
"""
CapG91# = Variant()
CapG92# = Variant()
CapG93# = Variant()
Cour# = Variant()
DRBar1# = Variant()
DRBar2# = Variant()
DRBar3# = Variant()
DTau# = Variant()
Speed1 = Variant()
jMax1% = Variant()
jMax2% = Variant()
jMax3% = Variant()
iLoop% = Variant()
iPrint% = Variant()
Beta21# = Variant()
Phi21# = Variant()
Bata31# = Variant()
Phi31# = Variant()
RA# = Variant()
RB# = Variant()
RC# = Variant()
RLR# = Variant()
RLT# = Variant()
Tau# = Variant()
TauMax# = Variant()
TPrint# = Variant()
U1# = vbObjectInitialize((3, 2798,), Variant)
UP1# = vbObjectInitialize((3, 2798,), Variant)
U2# = vbObjectInitialize((3, 2134,), Variant)
UP2# = vbObjectInitialize((3, 2134,), Variant)
U3# = vbObjectInitialize((3, 16543,), Variant)
UP3# = vbObjectInitialize((3, 16543,), Variant)
def CCyl1():
Pois1# = Variant()
Pois2# = Variant()
Pois3# = Variant()
Beta21# = Variant()
Beta31# = Variant()
SigRA# = Variant()
SigRD# = Variant()
Speed1# = Variant()
UOld11# = Variant()
Eigen1# = Variant()
Eigen2# = Variant()
Eigen3# = Variant()
EigMx1# = Variant()
EigMx2# = Variant()
EigMx3# = Variant()
RLTSq# = Variant()
EigSq# = Variant()
DTau1# = Variant()
DTau2# = Variant()
DTau3# = Variant()
D1# = Variant()
D2# = Variant()
D3# = Variant()
PartA# = Variant()
PartB# = Variant()
Hold1# = Variant()
Hold2# = Variant()
Hold3# = Variant()
VAvg# = Variant()
DLmdaT# = Variant()
RBar# = Variant()
RCauch# = Variant()
j% = Variant()
k% = Variant()
# ---------------------------------------------------------------------
# Converted from CCYL1.FOR (in P:\Archive\PHD1988\CCyl1993\CCYL1&2)
# Concentric Cylinders In An Unbounded Medium
# RA=1 RB=1.5
# 1 -------------- jMax1 RC=2 RD
# 1 -------------- jMax2 INFINITE
# 1 --------------------- jMax3
# PT Janele 14-10-04
# ---------------------------------------------------------------------
#COMMON /COM1/ CapG91, CapG92, CapG93
#COMMON /COM2/ Cour, DRBAR1, DRBAR2, DRBAR3, DTAU, SPEED1
#COMMON /COM3/ jMax1, jMax2, jMax3, ILOOP, IPRINT
#COMMON /COM4/ BETA21, PHI21, BETA31, PHI31
#COMMON /COM5/ RA, RB, RC, RLR, RLT
#COMMON /COM6/ TAU, TAUMAX, TPRINT
#COMMON /COM7/ U1(3,2798), UP1(3,2798)
#COMMON /COM8/ U2(3,2134), UP2(3,2134)
#COMMON /COM9/ U3(3,16543), UP3(3,16543)
#... OPEN OUTPUT DATA FILES
#OPEN (UNIT=2, FILE='CCYL1-2.DAT', STATUS='UNKNOWN')
#OPEN (UNIT=3, FILE='CCYL1-3.DAT', STATUS='UNKNOWN')
#OPEN (UNIT=4, FILE='CCYL1-4.DAT', STATUS='UNKNOWN')
#OPEN (UNIT=7, FILE='CCYL1-7.DAT', STATUS='UNKNOWN')
#OPEN (UNIT=8, FILE='CCYL1-8.DAT', STATUS='UNKNOWN')
#OPEN (UNIT=9, FILE='CCYL1-9.DAT', STATUS='UNKNOWN')
VBFiles.openFile(2, 'CCyl1-2Out.dat', 'w')
VBFiles.openFile(3, 'CCyl1-3Out.dat', 'w')
VBFiles.openFile(4, 'CCyl1-4Out.dat', 'w')
VBFiles.openFile(7, 'CCyl1-7Out.dat', 'w')
VBFiles.openFile(8, 'CCyl1-8Out.dat', 'w')
VBFiles.openFile(9, 'CCyl1-9Out.dat', 'a')
#... Define Data For Dimensionless Formulation
Cour = 0.999
DRBar3 = 0.0002
iPrint = 0
Pois1 = 0.3
Pois2 = 0.4
Pois3 = 0.35
Beta21 = 1.2
Phi21 = 1.2
Beta31 = 1.1
Phi31 = 1.1
RA = 1
RB = 1.5
RC = 2
SigRA = - 1
SigRD = 0
#... Initialize The Constants CapG91#, CapG92#, CapG93#
CapG91 = 2 * ( 1 + Pois1 ) / ( 3 * ( 1 - 2 * Pois1 ) ) / 9
CapG92 = 2 * ( 1 + Pois2 ) / ( 3 * ( 1 - 2 * Pois2 ) ) / 9
CapG93 = 2 * ( 1 + Pois3 ) / ( 3 * ( 1 - 2 * Pois3 ) ) / 9
#... Initialize Time Parameters
Tau = 0
Speed1 = Sqr(( 1 / 3 ) + 1 + CapG91 * 9)
TPrint = 0.125
TauMax = 0.75
#... Define Initial Eigenvalues
Eigen1 = Sqr(( 1 / 3 ) + 1 + CapG91 * 9)
Eigen2 = Sqr(( Phi21 / Beta21 ) * ( ( 1 / 3 ) + 1 + CapG92 * 9 ))
Eigen3 = Sqr(( Phi31 / Beta31 ) * ( ( 1 / 3 ) + 1 + CapG93 * 9 ))
Debug.Print(Eigen1, Eigen2, Eigen3)
# WRITE(6,100) Eigen1#, Eigen2#, Eigen3#
# 100 FORMAT(' ','INITIAL : Eigen1# Eigen2# Eigen3# = ',3(1X,F16.6))
#... Define DRBar1#, jMax1% AND jMax2% : (Material #1)
DRBar1 = DRBar3 * ( Eigen1 / Eigen3 )
DRBar2 = DRBar3 * ( Eigen2 / Eigen3 )
jMax1 = ( RB - RA ) / DRBar1 + 1.05
jMax2 = ( RC - RB ) / DRBar2 + 1.05
jMax3 = 16543
Debug.Print(jMax1, jMax2, jMax3)
# WRITE(6,120) jMax1%, jMax2%, jMax3%
# 120 FORMAT(' ','jMax1% jMax2% jMax3% = ',3I9)
#... Modify DRBar1#, DRBar2# So That(jMax1%-1)*DRBar = (RB#-RA#) ...
DRBar1 = ( RB - RA ) / ( jMax1 - 1 )
DRBar2 = ( RC - RB ) / ( jMax2 - 1 )
Debug.Print(DRBar1, DRBar2, DRBar3)
# WRITE(6,150) DRBar1#, DRBar2#, DRBar3#
# 150 FORMAT(' ','DRBar1# DRBar2# DRBar3# = ',3(1X,F16.6),//)
#... Define Initial Time Step Size : Use 10. As "Safety" Factor
DTau1 = Cour * DRBar1 / Eigen1
DTau2 = Cour * DRBar2 / Eigen2
DTau3 = Cour * DRBar3 / Eigen3
DTau = Application.WorksheetFunction.Min(DTau1, DTau2, DTau3) / 10
#... Material Is Initially In An Unstressed State
for j% in vbForRange(1, jMax1):
U1[1, j] = 0
U1[2, j] = 1
U1[3, j] = 1
for j% in vbForRange(1, jMax2):
U2[1, j] = 0
U2[2, j] = 1
U2[3, j] = 1
for j% in vbForRange(1, jMax3):
U3[1, j] = 0
U3[2, j] = 1
U3[3, j] = 1
# ======================================================================
#... Top Of Outer Time Step Loop : Do 800 Loop Here
for iLoop% in vbForRange(2, 30000):
#... Loop On Outer Surrounding Material Dependent On Wave Front; iloop%
jMax3 = iLoop + 2
if ( jMax3 > 16543 ) :
break
Tau = Tau + DTau
if ( Tau > TauMax ) :
break
Debug.Print(iLoop, Tau, jMax3)
#... Find New Lamda-R At Inside Boundary
RLT = U1(3, 1)
Call(FndRLR(SigRA, CapG91, 1))
U1[2, 1] = RLR
#... Search For Largest Eigenvalue And Define Time Step
EigMx1 = 0
for j% in vbForRange(1, jMax1):
RLR = U1(2, j)
RLT = U1(3, j)
RLTSq = RLT * RLT
EigSq = 1 + ( RLTSq * ( RLR * RLT ) ** ( - 4 / 3 ) ) / 3
EigSq = EigSq + CapG91 * ( 10 * RLR ** ( - 11 ) * RLT ** ( - 9 ) - RLR ** ( - 2 ) )
if ( EigSq >= EigMx1 ) :
EigMx1 = EigSq
EigMx2 = 0
for j% in vbForRange(1, jMax2):
RLR = U2(2, j)
RLT = U2(3, j)
RLTSq = RLT * RLT
EigSq = 1 + ( RLTSq * ( RLR * RLT ) ** ( - 4 / 3 ) ) / 3
EigSq = EigSq + CapG92 * ( 10 * RLR ** ( - 11 ) * RLT ** ( - 9 ) - RLR ** ( - 2 ) )
EigSq = ( Phi21 / Beta21 ) * EigSq
if ( EigSq >= EigMx2 ) :
EigMx2 = EigSq
EigMx3 = 0
for j% in vbForRange(1, jMax3):
RLR = U3(2, j)
RLT = U3(3, j)
RLTSq = RLT * RLT
EigSq = 1 + ( RLTSq * ( RLR * RLT ) ** ( - 4 / 3 ) ) / 3
EigSq = EigSq + CapG93 * ( 10 * RLR ** ( - 11 ) * RLT ** ( - 9 ) - RLR ** ( - 2 ) )
EigSq = ( Phi31 / Beta31 ) * EigSq
if ( EigSq >= EigMx3 ) :
EigMx3 = EigSq
DTau1 = Cour * DRBar1 / Sqr(EigMx1)
DTau2 = Cour * DRBar2 / Sqr(EigMx2)
DTau3 = Cour * DRBar3 / Sqr(EigMx3)
DTau = Application.WorksheetFunction.Min(DTau1, DTau2, DTau3)
D1 = DTau / DRBar1
D2 = DTau / DRBar2
D3 = DTau / DRBar3
# --------------------- Mac Cormack Predictor ------------------------
#... Special Case At j% = 1 : Material #1
RBar = RA
RLR = U1(2, 2)
RLT = U1(3, 2)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = PartA + CapG91 * PartB
RLR = U1(2, 1)
RLT = U1(3, 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = PartA + CapG91 * PartB
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( PartA + CapG91 * PartB ) / RBar
UP1[1, 1] = U1(1, 1) - D1 * ( Hold1 - Hold2 ) - DTau * Hold3
VAvg = 0.5 * ( U1(1, 1) + UP1(1, 1) )
DLmdaT = VAvg * DTau / RBar
UP1[3, 1] = U1(3, 1) + DLmdaT
RLT = UP1(3, 1)
Call(FndRLR(SigRA, CapG91, 1))
UP1[2, 1] = RLR
#... Predictor Step For 2 < j% < jMax1%-1 : Material #1
for j% in vbForRange(2, jMax1 - 1):
RBar = RA + ( j - 1 ) * DRBar1
RLR = U1(2, j + 1)
RLT = U1(3, j + 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = PartA + CapG91 * PartB
RLR = U1(2, j)
RLT = U1(3, j)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = PartA + CapG91 * PartB
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( PartA + CapG91 * PartB ) / RBar
UP1[1, j] = U1(1, j) - D1 * ( Hold1 - Hold2 ) - DTau * Hold3
UP1[2, j] = U1(2, j) + D1 * ( U1(1, j + 1) - U1(1, j) )
UP1[3, j] = U1(3, j) + DTau * U1(1, j) / RBar
#... Predictor : Special Case At j% = jMax1% : Material #1
# Forward Diff Scheme; Need Values In Material #2 At Grid Node 2
# Find Cauchy Radial Stress In Material #2; Grid Node 2
RLR = U2(2, 2)
RLT = U2(3, 2)
RCauch = RLR - RLT * ( RLR * RLT ) ** ( - 1 / 3 )
RCauch = RCauch + CapG92 * ( 1 / RLR - RLR ** ( - 10 ) * RLT ** ( - 9 ) )
RCauch = Phi21 * ( RCauch / RLT )
# Find Psuedo Rlr# Which Would Be In Material #1 And At Above Grid Node
# Retain Psuedo Rlt# In Material #1 = Actual Rlt# In Material #2
#Note That Resultant Rlr# Is Returned In Common Block
Call(FndRLR(RCauch, CapG91, 1))
# continue With Foreward Difference Scheme; Required RLr# Was Returned
RBar = RB
RLT = U2(3, 2)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = PartA + CapG91 * PartB
RLR = U1(2, jMax1)
RLT = U1(3, jMax1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = PartA + CapG91 * PartB
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( PartA + CapG92 * PartB ) / RBar
UP1[1, jMax1] = U1(1, jMax1) - D1 * ( Hold1 - Hold2 ) - DTau * Hold3
UP1[2, jMax1] = U1(2, jMax1) + D1 * ( U2(1, 2) - U1(1, jMax1) )
UP1[3, jMax1] = U1(3, jMax1) + DTau * U1(1, jMax1) / RBar
#... PREDICTOR STEP FOR 1 < j% < jMax2%-1 : MATERIAL #2
for j% in vbForRange(1, jMax2 - 1):
RBar = RB + ( j - 1 ) * DRBar2
RLR = U2(2, j + 1)
RLT = U2(3, j + 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
RLR = U2(2, j)
RLT = U2(3, j)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB ) / RBar
UP2[1, j] = U2(1, j) - D2 * ( Hold1 - Hold2 ) - DTau * Hold3
UP2[2, j] = U2(2, j) + D2 * ( U2(1, j + 1) - U2(1, j) )
UP2[3, j] = U2(3, j) + DTau * U2(1, j) / RBar
#... Predictor : Special Case At j% = jMax2% : Material #2
# Forward Diff Scheme; Need Values In Material #3 At Grid Node 2
# Find Cauchy Ra#Dial Stress In Material #3; Grid Node 2
RLR = U3(2, 2)
RLT = U3(3, 2)
RCauch = RLR - RLT * ( RLR * RLT ) ** ( - 1 / 3 )
RCauch = RCauch + CapG93 * ( 1 / RLR - RLR ** ( - 10 ) * RLT ** ( - 9 ) )
RCauch = Phi31 * ( RCauch / RLT )
# Find Psuedo Rlr# Which Would Be In Material #2 And At Above Grid Node
# Retain Psuedo Rlt# In Material #1 = Actual Rlt# In Material #2
# Note That Resultant Rlr# Is Returned In Common Block
Call(FndRLR(RCauch, CapG92, Phi21))
# Continue With Foreward Difference Scheme; Required RLr# Was Returned
RBar = RC
RLT = U3(3, 2)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
RLR = U2(2, jMax2)
RLT = U2(3, jMax2)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB ) / RBar
UP2[1, jMax2] = U2(1, jMax2) - D2 * ( Hold1 - Hold2 ) - DTau * Hold3
UP2[2, jMax2] = U2(2, jMax2) + D2 * ( U3(1, 2) - U2(1, jMax2) )
UP2[3, jMax2] = U2(3, jMax2) + DTau * U2(1, jMax2) / RBar
#... Continuation Of Predictor Step For 1 < j% < jMax3%-1 : Material #3
for j% in vbForRange(1, jMax3 - 1):
RBar = RC + ( j - 1 ) * DRBar3
RLR = U3(2, j + 1)
RLT = U3(3, j + 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB )
RLR = U3(2, j)
RLT = U3(3, j)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB )
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB ) / RBar
UP3[1, j] = U3(1, j) - D3 * ( Hold1 - Hold2 ) - DTau * Hold3
UP3[2, j] = U3(2, j) + D3 * ( U3(1, j + 1) - U3(1, j) )
UP3[3, j] = U3(3, j) + DTau * U3(1, j) / RBar
# --------------------- MACORMACK CORRECTOR ------------------------
#... Corrector : Special Case At J% = 1 : Material #1
#... Use Foreward Difference Here
RBar = RA
RLR = UP1(2, 2)
RLT = UP1(3, 2)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = PartA + CapG91 * PartB
RLR = UP1(2, 1)
RLT = UP1(3, 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = PartA + CapG91 * PartB
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( PartA + CapG91 * PartB ) / RBar
UOld11 = U1(1, 1)
U1[1, 1] = 0.5 * ( UP1(1, 1) + U1(1, 1) - D1 * ( Hold1 - Hold2 ) - DTau * Hold3 )
VAvg = 0.5 * ( U1(1, 1) + UOld11 )
DLmdaT = VAvg * DTau / RBar
U1[3, 1] = U1(3, 1) + DLmdaT
RLT = U1(3, 1)
Call(FndRLR(SigRA, CapG91, 1))
U1[2, 1] = RLR
#... Corrector Step For 2 < j% < jMax1% : Material #1
for j% in vbForRange(2, jMax1):
RBar = RA + ( j - 1 ) * DRBar1
RLR = UP1(2, j)
RLT = UP1(3, j)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = PartA + CapG91 * PartB
RLR = UP1(2, j - 1)
RLT = UP1(3, j - 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = PartA + CapG91 * PartB
RLR = UP1(2, j)
RLT = UP1(3, j)
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( PartA + CapG91 * PartB ) / RBar
U1[1, j] = 0.5 * ( UP1(1, j) + U1(1, j) - D1 * ( Hold1 - Hold2 ) - DTau * Hold3 )
U1[2, j] = 0.5 * ( UP1(2, j) + U1(2, j) + D1 * ( UP1(1, j) - UP1(1, j - 1) ) )
U1[3, j] = 0.5 * ( UP1(3, j) + U1(3, j) + DTau * UP1(1, j) / RBar )
#... Corrector : Special Case At jMax2% = 1 : Material #2
# Backward Diff Scheme; Need Values In Material #1 At Grid Node jMax1%-1
# Find Cauchy Radial Stress In Material #1; Grid Node jMax1%-1
RLR = UP1(2, jMax1 - 1)
RLT = UP1(3, jMax1 - 1)
RCauch = RLR - RLT * ( RLR * RLT ) ** ( - 1 / 3 )
RCauch = RCauch + CapG91 * ( 1 / RLR - RLR ** ( - 10 ) * RLT ** ( - 9 ) )
RCauch = RCauch / RLT
# Find Psuedo Rlr# Which Would Be In Material #1 And At Above Grid Node
# Retain Psuedo Rlt# In Material #1 = Actual Rlt# In Material #2
# Note That Resultant Rlr# Is Returned In Common Block
Call(FndRLR(RCauch, CapG92, Phi21))
# Continue With Backward Difference Scheme; Required RLr# Was Returned
RBar = RB
RLT = UP1(3, jMax1 - 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
RLR = UP2(2, 1)
RLT = UP2(3, 1)
# ptj fails here and below... but problem likely to be in Predictor step
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB ) / RBar
U2[1, 1] = 0.5 * ( UP2(1, 1) + U2(1, 1) - D1 * ( Hold1 - Hold2 ) - DTau * Hold3 )
U2[2, 1] = UP2(2, 1) + U2(2, 1) + D1 * ( UP2(1, 1) - UP1(1, jMax1 - 1) )
U2[2, 1] = 0.5 * U2(2, 1)
U2[3, 1] = 0.5 * ( UP2(3, 1) + U2(3, 1) + DTau * UP2(1, 1) / RBar )
#... Continuation Of Corrector Step For 2 < j% < jMax2% : Material #2
for j% in vbForRange(2, jMax2):
RBar = RB + ( j - 1 ) * DRBar2
RLR = UP2(2, j)
RLT = UP2(3, j)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
RLR = UP2(2, j - 1)
RLT = UP2(3, j - 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB )
RLR = UP2(2, j)
RLT = UP2(3, j)
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi21 / Beta21 ) * ( PartA + CapG92 * PartB ) / RBar
U2[1, j] = 0.5 * ( UP2(1, j) + U2(1, j) - D2 * ( Hold1 - Hold2 ) - DTau * Hold3 )
U2[2, j] = 0.5 * ( UP2(2, j) + U2(2, j) + D2 * ( UP2(1, j) - UP2(1, j - 1) ) )
U2[3, j] = 0.5 * ( UP2(3, j) + U2(3, j) + DTau * UP2(1, j) / RBar )
#... Corrector : Special Case At Jmax3% = 1 : Material #3
# Backward Diff Scheme; Need Values In Material #2 At Grid Node Jmax2%-2
# Find Cauchy Ra#Dial Stress In Material #2; Grid Node Jmax2%-1
RLR = UP2(2, jMax2 - 1)
RLT = UP2(3, jMax2 - 1)
RCauch = RLR - RLT * ( RLR * RLT ) ** ( - 1 / 3 )
RCauch = RCauch + CapG92 * ( 1 / RLR - RLR ** ( - 10 ) * RLT ** ( - 9 ) )
RCauch = Phi21 * ( RCauch / RLT )
# Find Psuedo Rlr# Which Would Be In Material #1 And At Above Grid Node
# Retain Psuedo Rlt# In Material #1 = Actual Rlt# In Material #2
# Note That Resultant Rlr# Is Returned In Common Block
Call(FndRLR(RCauch, CapG93, Phi31))
# Continue With Backward Difference Scheme; Required Rlr# Was Returned
RBar = RC
RLT = UP2(3, jMax2 - 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB )
RLR = UP3(2, 1)
RLT = UP3(3, 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB )
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB ) / RBar
U3[1, 1] = 0.5 * ( UP3(1, 1) + U3(1, 1) - D2 * ( Hold1 - Hold2 ) - DTau * Hold3 )
U3[2, 1] = UP3(2, 1) + U3(2, 1) + D2 * ( UP3(1, 1) - UP2(1, jMax2 - 1) )
U3[2, 1] = 0.5 * U3(2, 1)
U3[3, 1] = 0.5 * ( UP3(3, 1) + U3(3, 1) + DTau * UP3(1, 1) / RBar )
#... Continuation Of Corrector Step For 2 < J% < Jmax3%-1 : Material #3
for j% in vbForRange(2, jMax3 - 1):
RBar = RC + ( j - 1 ) * DRBar3
RLR = UP3(2, j)
RLT = UP3(3, j)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold1 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB )
RLR = UP3(2, j - 1)
RLT = UP3(3, j - 1)
PartA = RLT * ( RLR * RLT ) ** ( - 1 / 3 ) - RLR
PartB = RLT ** ( - 9 ) * RLR ** ( - 10 ) - 1 / RLR
Hold2 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB )
RLR = UP3(2, j)
RLT = UP3(3, j)
PartA = ( RLT - RLR ) * ( 1 + ( RLR * RLT ) ** ( - 1 / 3 ) )
PartB = ( 1 / RLT - 1 / RLR ) * ( 1 - ( RLR * RLT ) ** ( - 9 ) )
Hold3 = ( Phi31 / Beta31 ) * ( PartA + CapG93 * PartB ) / RBar
U3[1, j] = 0.5 * ( UP3(1, j) + U3(1, j) - D3 * ( Hold1 - Hold2 ) - DTau * Hold3 )
U3[2, j] = 0.5 * ( UP3(2, j) + U3(2, j) + D3 * ( UP3(1, j) - UP3(1, j - 1) ) )
U3[3, j] = 0.5 * ( UP3(3, j) + U3(3, j) + DTau * UP3(1, j) / RBar )
# ----------------- Write Out Results And Loop Time -------------------
#... Print Out Every 10Th Data Point
# WRITE(7,750) iLoop%, Tau#, Tau#*Speed1#, U1#(1,1), U1#(2,1), U1#(3,1)
VBFiles.writeText(7, iLoop, "\t", Tau, "\t", Tau * Speed1, "\t", U1(1, 1), "\t", U1(2, 1), "\t", U1(3, 1), '\n')
# WRITE(2,750) iLoop%, Tau#, Tau#*Speed1#, U1#(1,jMax1%), U1#(2,jMax1%), U1#(3,jMax1%)
VBFiles.writeText(2, iLoop, "\t", Tau, "\t", Tau * Speed1, "\t", U1(1, jMax1), "\t", U1(2, jMax1), "\t", U1(3, jMax1), '\n')
# WRITE(3,750) iLoop%, Tau#, Tau#*Speed1#, U2#(1,1), U2#(2,1), U2#(3,1)
VBFiles.writeText(3, iLoop, "\t", Tau, "\t", Tau * Speed1, "\t", U2(1, 1), "\t", U2(2, 1), "\t", U2(3, 1), '\n')
# WRITE(8,750) iLoop%, Tau#, Tau#*Speed1#, U2#(1,jMax2%), U2#(2,jMax2%), U2#(3,jMax2%)
VBFiles.writeText(8, iLoop, "\t", Tau, "\t", Tau * Speed1, "\t", U2(1, jMax2), "\t", U2(2, jMax2), "\t", U2(3, jMax2), '\n')
# WRITE(4,750) iLoop%, Tau#, Tau#*Speed1#, U3#(1,1), U3#(2,1), U3#(3,1)
VBFiles.writeText(4, iLoop, "\t", Tau, "\t", Tau * Speed1, "\t", U3(1, 1), "\t", U3(2, 1), "\t", U3(3, 1), '\n')
k = Int(( Tau + DTau * ( 1 + TPrint ) ) / TPrint)
if ( ( Tau <= k * TPrint ) and ( Tau + DTau ) > k * TPrint ) :
Call(FShock)
#... Bottom Of Outer Time Step Loop
# close all the output files
VBFiles.closeFile(2)
VBFiles.closeFile(3)
VBFiles.closeFile(4)
VBFiles.closeFile(7)
VBFiles.closeFile(8)
VBFiles.closeFile(9)
def FndRLR(Target#, CAPG9X#, PhiXX#):
i% = Variant()
RLTSq# = Variant()
RCauch# = Variant()
Slope# = Variant()
ConvTol# = Variant()
# COMMON /COM5/ RA, RB, RC, RLR#, RLT#
#... Start With Guessed Value Of RLR# and Loop for Solution
RLR = 0.8
RLTSq = RLT * RLT
#ConvTol# = 10 ^ (-6)
ConvTol = 10 ** ( - 12 )
for i% in vbForRange(1, 2000):
# converge via Newton method
RCauch = RLR - RLT * ( RLR * RLT ) ** ( - 1 / 3 )
RCauch = RCauch + CAPG9X * ( 1 / RLR - RLR ** ( - 10 ) * RLT ** ( - 9 ) )
RCauch = PhiXX * ( RCauch / RLT )
if ( Abs(RCauch - Target) < ConvTol ) :
break
Slope = 1 + RLTSq * ( ( RLR * RLT ) ** ( - 4 / 3 ) ) / 3
Slope = Slope + CAPG9X * ( 10 * RLR ** ( - 11 ) * RLT ** ( - 9 ) - RLR ** ( - 2 ) )
Slope = PhiXX * ( Slope / RLT )
RLR = RLR - ( RCauch - Target ) / Slope
# Debug.Print i%, RLR#, RCauch# - Target#
#WRITE(6,1020) k% - 1
#1020 FORMAT(' ','ERROR IN FNDRLR# AFTER ',I8,' ITERATIONS')
#Stop
def FShock():
k% = Variant()
RAvg# = Variant()
Radius# = Variant()
YAvg# = Variant()
Total# = Variant()
# COMMON /COM1/ CAPG91, CAPG92, CAPG93
# COMMON /COM2/ COUR, DRBar1#, DRBAR2, DRBAR3, DTau#, SPEED1
# COMMON /COM3/ JMAX1, JMAX2, JMAX3, ILOOP, IPRINT
# COMMON /COM5/ RA, RB, RC, RLR, RLT
# COMMON /COM6/ Tau, TAUMAX, TPRINT
# COMMON /COM7/ U1(3,2798), UP1(3,2798)
# COMMON /COM8/ U3(3,2134), UP2(3,2134)
# COMMON /COM9/ U3(3,16543), UP3(3,16543)
# OPEN (UNIT=9, FILE='CCYL1-9.DAT', STATUS='UNKNOWN')
#Open "CCyl1-9Out.dat" For Append As #9
VBFiles.writeText(9, '-----------', '\n')
VBFiles.writeText(9, 'Solution Matrix U(i%,k%) At Time Tau# = ', Tau, '\n')
# Write(9,1200)
# 1200 FORMAT(///)
# Write(9,1220) Tau#
# 1220 FORMAT(/,'Solution Matrix U1#(I,k%) At Time Tau# = ',E16.9,/)
Total = 0
for k% in vbForRange(1, jMax1):
RAvg = RA + ( k - 1 ) * DRBar1 + DRBar1 / 2
Radius = RA + ( k - 1 ) * DRBar1
YAvg = 0.5 * ( U1(1, k) + U1(1, k + 1) )
Total = Total + YAvg * RAvg * RAvg * DRBar1
VBFiles.writeText(9, k, "\t", Radius, "\t", U1(1, k), "\t", U1(2, k), "\t", U1(3, k), '\n')
#WRITE(9,1240) k%, Radius#, U1#(1,k%), U1#(2,k%), U1#(3,k%)
#1240 FORMAT(' ',I9,4(1X,E16.9))
for k% in vbForRange(1, jMax2):
RAvg = RB + ( k - 1 ) * DRBar2 + DRBar2 / 2
Radius = RB + ( k - 1 ) * DRBar2
YAvg = 0.5 * ( U2(1, k) + U2(1, k + 1) )
Total = Total + YAvg * RAvg * RAvg * DRBar2
VBFiles.writeText(9, k, "\t", Radius, "\t", U2(1, k), "\t", U2(2, k), "\t", U2(3, k), '\n')
for k% in vbForRange(1, jMax3):
RAvg = RC + ( k - 1 ) * DRBar3 + DRBar3 / 2
Radius = RC + ( k - 1 ) * DRBar3
YAvg = 0.5 * ( U3(1, k) + U3(1, k + 1) )
Total = Total + YAvg * RAvg * RAvg * DRBar3
VBFiles.writeText(9, k, "\t", Radius, "\t", U3(1, k), "\t", U3(2, k), "\t", U3(3, k), '\n')
#WRITE(9,1300) JMAX1, JMAX2, JMAX3, TAU, TOTAL
# 1300 FORMAT(' ','JMAX1 JMAX2 JMAX3 TAU TOTAL ',3I7,2(1X,E16.9))
VBFiles.writeText(9, jMax1, "\t", jMax2, "\t", jMax3, "\t", Tau, "\t", Total, '\n')
#Close #9
# VB2PY (UntranslatedCode) Attribute VB_Name = "CCylMain"
# VB2PY (UntranslatedCode) Option Explicit
# VB2PY (UntranslatedCode) Option Base 1
Write, Run & Share Python code online using OneCompiler's Python online compiler for free. It's one of the robust, feature-rich online compilers for python language, supporting both the versions which are Python 3 and Python 2.7. Getting started with the OneCompiler's Python editor is easy and fast. The editor shows sample boilerplate code when you choose language as Python or Python2 and start coding.
OneCompiler's python online editor supports stdin and users can give inputs to programs using the STDIN textbox under the I/O tab. Following is a sample python program which takes name as input and print your name with hello.
import sys
name = sys.stdin.readline()
print("Hello "+ name)
Python is a very popular general-purpose programming language which was created by Guido van Rossum, and released in 1991. It is very popular for web development and you can build almost anything like mobile apps, web apps, tools, data analytics, machine learning etc. It is designed to be simple and easy like english language. It's is highly productive and efficient making it a very popular language.
When ever you want to perform a set of operations based on a condition IF-ELSE is used.
if conditional-expression
#code
elif conditional-expression
#code
else:
#code
Indentation is very important in Python, make sure the indentation is followed correctly
For loop is used to iterate over arrays(list, tuple, set, dictionary) or strings.
mylist=("Iphone","Pixel","Samsung")
for i in mylist:
print(i)
While is also used to iterate a set of statements based on a condition. Usually while is preferred when number of iterations are not known in advance.
while condition
#code
There are four types of collections in Python.
List is a collection which is ordered and can be changed. Lists are specified in square brackets.
mylist=["iPhone","Pixel","Samsung"]
print(mylist)
Tuple is a collection which is ordered and can not be changed. Tuples are specified in round brackets.
myTuple=("iPhone","Pixel","Samsung")
print(myTuple)
Below throws an error if you assign another value to tuple again.
myTuple=("iPhone","Pixel","Samsung")
print(myTuple)
myTuple[1]="onePlus"
print(myTuple)
Set is a collection which is unordered and unindexed. Sets are specified in curly brackets.
myset = {"iPhone","Pixel","Samsung"}
print(myset)
Dictionary is a collection of key value pairs which is unordered, can be changed, and indexed. They are written in curly brackets with key - value pairs.
mydict = {
"brand" :"iPhone",
"model": "iPhone 11"
}
print(mydict)
Following are the libraries supported by OneCompiler's Python compiler
| Name | Description |
|---|---|
| NumPy | NumPy python library helps users to work on arrays with ease |
| SciPy | SciPy is a scientific computation library which depends on NumPy for convenient and fast N-dimensional array manipulation |
| SKLearn/Scikit-learn | Scikit-learn or Scikit-learn is the most useful library for machine learning in Python |
| Pandas | Pandas is the most efficient Python library for data manipulation and analysis |
| DOcplex | DOcplex is IBM Decision Optimization CPLEX Modeling for Python, is a library composed of Mathematical Programming Modeling and Constraint Programming Modeling |