|
#============================================================================== |
|
# THIS PROGRAM CALCULATES THE DEFORMATION OF A THIN ELASTIC PLATE |
|
# IN RESPONSE TO DISTRIBUTED LOADS. THE PLATE |
|
# THICKNESS, AND HENCE THE FLEXURAL RIGIDITY MAY VARY SPATIALLY. |
|
# THE PLATE MAY BE EITHER CONTINUOUS OR BROKEN. |
|
# THE LEFT SIDE (COL=3) IS CONSIDERED TO BE SIDE 1, |
|
# THE RIGHT SIDE (COL=NCOLS-2) IS SIDE 2, |
|
# THE BOTTOM (ROW=3) IS SIDE 3, AND |
|
# THE TOP (ROW=NROWS-2) IS SIDE 4. |
|
# SEE TIMOSHENKO & WOINOWSKY-KRIEGER "THEORY OF PLATES & SHELLS" |
|
# 1959, PG 174 FOR DifFERENTIAL EQUATION. |
|
# SOLUTION OF DifFERENTIAL EQUATION IS BY |
|
# SIMULTANEOUS OVER-RELAXATION GAUSS-SEIDEL ITERATION. |
|
# THE MAXIMUM ARRAY SIZE IS 100X100 GRIDS BUT A SUB-SET OF THIS ARRAY |
|
# MAY BE USED TO DECREASE EXECUTION TIME. |
|
# THE ENTIRE ARRAY AREA USED IS WRITTEN TO AN OUTPUT FILE FROM WHICH |
|
# SPECifIC CONTOUR AND CROSS-SECTION PLOTS MAY BE SELECTED AFTER THE |
|
# JOB IS COMPLETE. |
|
|
|
#============================================================================== |
|
# WRITTEN BY GARRY QUINLAN, FEB. 1990 |
|
|
|
# LAST MODifICATION - MAY 12, 1990 |
|
|
|
# PURPOSE OF MODifICATION - TO CORRECT (I HOPE#) THE EXPRESSIONS FOR |
|
# BOUNDARY CONDITIONS AT A BROKEN EDGE |
|
#============================================================================== |
|
|
|
#============================================================================== |
|
# # Converted to python by Michael Clark 2013 |
|
# Made slightly nicer for python by Michael Clark 2014 |
|
#============================================================================== |
|
|
|
import numpy # or numpypy for pypy |
|
|
|
|
|
|
|
def flex3db (D, P, GAMMA, dx, ibpsw, nits, tol, pr = 0.25, omega = 1.7, fx = 0, fy = 0, fxy = 0): |
|
# A THIN ELASTIC PLATE HAS SPATIALLY VARYING FLEXURAL RIGIDITY |
|
# SPECIfIED BY THE ARRAY D. EACH ELEMENT OF THIS ARRAY REPRESENTS |
|
# THE FLEXURAL RIGIDITY OF THE CORRESPONDING GRID WITHIN A 2-D |
|
# CARTESIAN COORDINATE SYSTEM. THE GRID SPACING IN THIS COORDINATE |
|
# SYSTEM IS THE 4TH ROOT OF THE VARIABLE DEL4 AND IS SPECifIED IN |
|
# METRES. |
|
|
|
# THE PLATE DESCRIBED BY THE FLEXURAL RIGIDITY ARRAY D IS SUBJECT |
|
# TO A SPATIALLY VARYING LOAD GIVEN BY THE ARRAY P. |
|
|
|
# THE DEFORMATION PRODUCED IN THE PLATE BY THE APPLIED LOAD IS |
|
# CALCULATED AND SAVED IN THE ARRAY W. |
|
|
|
# THE BOUYANCY FORCE PRODUCED BY THE DISPLACEMENT IS RELATED TO |
|
# THE DENSITY CONTRAST BETWEEN THE DISPLACED MANTLE AND WHATEVER |
|
# MATERIAL INFILLS THE DEPRESSION. THIS DENSITY CONTRAST IS |
|
# SPECifIED IN THE VARIABLE GAMMA = (DENSITY OF MANTLE - DENSITY |
|
# OF INFILL)*ACCELERATION DUE TO GRAVITY. |
|
|
|
# THE APPROPRIATE DifFERENTIAL EQUATION IS A 4TH ORDER PDE WHICH |
|
# IS DERIVED IN TIMOSHENKO & WOINOWSKY-KRIEGER "THEORY OF PLATES AND |
|
# SHELLS". |
|
|
|
# THE ASSUMED BOUNDARY CONDITIONS ARE ZERO DISPLACEMENT AT |
|
# THE EDGES OF THE GRID. THIS BOUNDARY VALUE PROBLEM IS |
|
# SOLVED BY SUCCESSIVE APPROXIMATIONS. THE def TERMINATES |
|
# WHEN EITHER: (1) THE MAXIMUM ALLOWABLE ITERATIONS (NITS) HAVE BEEN |
|
# MADE OR (2) WHEN THE CENTRAL ROW (IC) OF THE W ARRAY SATISFIES THE |
|
# GOVERNING PDE TO WITHIN A SPECifIED TOLERANCE (TOL). |
|
|
|
# IN SOME PROGRAMS USING THIS def THE FLEXURAL RIGIDITY WILL |
|
# BE THE SAME FROM ONE CALL OF THE def TO THE NEXT BUT IN |
|
# OTHERS IT MAY VARY. THE TERMS T1,T2,...T12 INVOLVE VARIOUS |
|
# COMBINATIONS OF PARTIAL DERIVATIVES OF THE FLEXURAL RIGIDITY. |
|
# THESE ARRAYS SHOULD BE CALCULATED ON ENTRY TO THE def ONLY |
|
# if FLEXURAL RIGIDITY HAS CHANGED SINCE THE LAST CALL TO THE |
|
# def. THIS NECESSITY IS IDENTifIED BY THE SWITCH ISW. |
|
# if ISW = 1, THE ARRAYS WILL BE CALCULATED; OTHERWISE WHATEVER |
|
# VALUES THEY CONTAIN WILL BE USED. |
|
|
|
# ****IT IS IMPORTANT TO NOTE THAT A PERIMETER 2 GRIDS WIDE MUST |
|
# BE INCLUDED IN THE SPECifICATION OF NROWS AND NCOLS TO ACCOMMODATE |
|
# THE BOUNDARY CONDITIONS FOR A CONTINUOUS PLATE I.E. THAT THE |
|
# DISPLACEMENT AND FIRST DERIVATIVE OF DISPLACEMENT VANISH. |
|
# THUS, if NROWS=20 AND NCOLS=15, THE TRUE WORKING AREA IS |
|
# ONLY 16X11.****** |
|
|
|
# THE PLATE MAY BE BROKEN ON 1 TO 4 SIDES. if BROKEN, THE BOUNDARY |
|
# CONDITIONS REQUIRE AN EXTRA COLUMN TO SPECifY THAT THE 2ND AND 3RD |
|
# DERIVATIVES VANISH. THIS FURTHER REDUCES THE SIZE OF THE TRUE WORKING |
|
# AREA. |
|
nrows, ncols = D.shape |
|
print "Running f3Dflex for {nits} iterations, and a {nrows}x{ncols} grid".format(nits = nits, nrows = nrows, ncols = ncols) |
|
converged = False |
|
|
|
|
|
T1 = numpy.zeros((nrows, ncols)) |
|
T2 = numpy.zeros((nrows, ncols)) |
|
T3 = numpy.zeros((nrows, ncols)) |
|
T4 = numpy.zeros((nrows, ncols)) |
|
T5 = numpy.zeros((nrows, ncols)) |
|
T6 = numpy.zeros((nrows, ncols)) |
|
T7 = numpy.zeros((nrows, ncols)) |
|
T8 = numpy.zeros((nrows, ncols)) |
|
T9 = numpy.zeros((nrows, ncols)) |
|
T10 = numpy.zeros((nrows, ncols)) |
|
T11 = numpy.zeros((nrows, ncols)) |
|
T12 = numpy.zeros((nrows, ncols)) |
|
T13 = numpy.zeros((nrows, ncols)) |
|
W = numpy.zeros((nrows, ncols)) |
|
|
|
onemom = 1.0 - omega |
|
twompr = 2.0 - pr |
|
onempr = 1.0 - pr |
|
prvint = 0.0 |
|
del4 = dx ** 4 |
|
delsq = numpy.sqrt(del4) |
|
|
|
# MULTIPLY FX,FY AND FXY BY DELSQ SO THAT ALL TERMS INVOLVING |
|
# PARTIAL DERIVATIVES HAVE A COMMON DENOMINATOR OF DEL4. |
|
fx = fx * delsq |
|
fy = fy * delsq |
|
fxy = fxy * delsq |
|
|
|
|
|
# DETERMINE THE TERMS INVOLVING PARTIAL DERIVATIVES OF FLEXURAL |
|
# RIGIDITY |
|
# NOTE THAT IN-PLANE FORCES MODifY TERMS T2,T4,T9,T11 AND T13 |
|
for i in xrange(3 - 1, nrows - 2): |
|
for j in xrange(3 - 1, ncols - 2): |
|
T1[i, j] = D[i + 1, j] |
|
T2[i, j] = D[i + 1, j] + D[i, j - 1] - (onempr / 8.0) * (D[i + 1, j + 1] - D[i - 1, j + 1] - D[i + 1, j - 1] + D[i - 1, j - 1]) - 2.0 * fxy |
|
T3[i, j] = (2.0 * onempr - 4.0) * D[i, j] - 4.0 * D[i + 1, j] - onempr * (D[i, j + 1] + D[i, j - 1]) - fx |
|
T4[i, j] = D[i + 1, j] + D[i, j + 1] + (onempr / 8.0) * (D[i + 1, j + 1] + D[i - 1, j - 1] - D[i - 1, j + 1] - D[i + 1, j - 1]) + 2.0 * fxy |
|
T5[i, j] = D[i, j - 1] |
|
T6[i, j] = (2.0 * onempr - 4.0) * D[i, j] - 4.0 * D[i, j - 1] - onempr * (D[i + 1, j] + D[i - 1, j]) - fy |
|
T7[i, j] = (2.0 * onempr - 4.0) * D[i, j] - 4.0 * D[i, j + 1] - onempr * (D[i + 1, j] + D[i - 1, j]) - fy |
|
T8[i, j] = D[i, j + 1] |
|
T9[i, j] = D[i - 1, j] + D[i, j - 1] + (onempr / 8.0) * (D[i + 1, j + 1] - D[i - 1, j + 1] - D[i + 1, j - 1] + D[i - 1, j - 1]) + 2.0 * fxy |
|
T10[i, j] = (2.0 * onempr - 4.0) * D[i, j] - 4.0 * D[i - 1, j] - onempr * (D[i, j + 1] + D[i, j - 1]) - fx |
|
T11[i, j] = D[i - 1, j] + D[i, j + 1] - (onempr / 8.0) * (D[i + 1, j + 1] - D[i - 1, j + 1] - D[i + 1, j - 1] + D[i - 1, j - 1]) - 2.0 * fxy |
|
T12[i, j] = D[i - 1, j] |
|
T13[i, j] = (1.0 + 2.0 * onempr) * (D[i + 1, j] + D[i - 1, j] + D[i, j + 1] + D[i, j - 1]) + (16.0 - 8.0 * onempr) * D[i, j] + del4 * GAMMA[i, j] + 2.0 * (fx + fy) |
|
|
|
|
|
|
|
for its in xrange(1, nits + 1): |
|
# ODD NUMBERED ITERATIONS LOOP THROUGH ARRAY FROM LOW TO HIGH |
|
# WHEREAS EVEN NUMBERED ITERATIONS LOOP FROM HIGH TO LOW. THIS |
|
# IS TO PREVENT THE SOLUTION FROM BECOMING MORE ACCURATE IN ONE |
|
# CORNER THAN THE OTHER AND SO, IN PRINCIPLE, TO SPEED CONVERGENCE. |
|
if (its % 2 == 0): |
|
for i in xrange(3 - 1, nrows - 2): |
|
for j in xrange(3 - 1, ncols - 2): |
|
residu = P[i, j] * del4 |
|
residu -= W[i + 2, j ] * T1[i, j] |
|
residu -= W[i + 1, j - 1] * T2[i, j] |
|
residu -= W[i + 1, j ] * T3[i, j] |
|
residu -= W[i + 1, j + 1] * T4[i, j] |
|
residu -= W[i , j - 2] * T5[i, j] |
|
residu -= W[i , j - 1] * T6[i, j] |
|
residu -= W[i , j + 1] * T7[i, j] |
|
residu -= W[i , j + 2] * T8[i, j] |
|
residu -= W[i - 1, j - 1] * T9[i, j] |
|
residu -= W[i - 1, j ] * T10[i, j] |
|
residu -= W[i - 1, j + 1] * T11[i, j] |
|
residu -= W[i - 2, j ] * T12[i, j] |
|
residu -= W[i , j ] * T13[i, j] |
|
W[i, j] += omega * residu / T13[i, j] |
|
else: |
|
for i in xrange(nrows - 3, 3 - 2, -1): |
|
for j in xrange(ncols - 3, 3 - 2, -1): |
|
residu = P[i, j] * del4 |
|
residu -= W[i + 2, j ] * T1[i, j] |
|
residu -= W[i + 1, j - 1] * T2[i, j] |
|
residu -= W[i + 1, j ] * T3[i, j] |
|
residu -= W[i + 1, j + 1] * T4[i, j] |
|
residu -= W[i , j - 2] * T5[i, j] |
|
residu -= W[i , j - 1] * T6[i, j] |
|
residu -= W[i , j + 1] * T7[i, j] |
|
residu -= W[i , j + 2] * T8[i, j] |
|
residu -= W[i - 1, j - 1] * T9[i, j] |
|
residu -= W[i - 1, j ] * T10[i, j] |
|
residu -= W[i - 1, j + 1] * T11[i, j] |
|
residu -= W[i - 2, j ] * T12[i, j] |
|
residu -= W[i , j ] * T13[i, j] |
|
W[i, j] += omega * residu / T13[i, j] |
|
|
|
|
|
|
|
|
|
# TEST CONVERGENCE OF SOLUTION. |
|
# CONVERGENCE CRITERION IS THAT THE INTEGRATED VALUE OF DISPLACEMENT |
|
# ALONG THE SPECifIED SECTION IS NO MORE THAN A SPECifIED AMOUNT |
|
# (TOL) DifFERENT FROM ITS VALUE ON PREVIOUS ITERATION (PRVINT). |
|
valint = numpy.sum(W) |
|
if valint != 0: |
|
change = numpy.abs((valint - prvint) / valint) |
|
else: |
|
change = 0 |
|
prvint = valint |
|
|
|
if (change < tol): |
|
print ' WITHIN TOLERANCE AFTER ', its, ' ITERATIONS', change |
|
converged = True |
|
break |
|
|
|
if its in xrange(0, nits, nits / 20): # 20 progress reports |
|
print " Iteration {its} of {nits}, change={change:f}".format(its = its, nits = nits, change = change) |
|
|
|
# if A SIDE OF THE STUDY AREA IS NOT BROKEN: |
|
# FORCE THE 1ST DERIVATIVES TO ZERO AT THE BOUNDARY. |
|
# THESE MODifICATIONS WILL AFFECT THE CALCULATIONS IN |
|
# THE NEXT ITERATION WHEN THESE MODifIED GRIDS ARE USED |
|
# TO DETERMINE DERIVATIVES. |
|
|
|
# if THE PROBLEM IS A BROKEN-PLATE ONE, THE BREAK IS ASSUMED TO BE |
|
# ALONG ONE OR MORE OF THE PLATE EDGES. |
|
# THE APPROPRIATE CONDITIONS ARE DISCUSSED IN TIMOSHENKO & W-K |
|
# PAGE 84. |
|
|
|
# NOTE THAT THE BROKEN PLATE PROBLEM REQUIRES AN EXTRA COLUMN |
|
# OUTBOARD OF THE AREA OF INTEREST TO ALLOW SETTING OF 3RD DERIVATIVES |
|
# ALONG THE BREAK. THIS IS ACCOMPLISHED BY TAKING COLUMN 3 RATHER THAN |
|
# COLUMN 2 TO BE THE START OF THE STUDY AREA. |
|
|
|
# SET BOUNDARY CONDITIONS FOR LEFT SIDE OF STUDY AREA |
|
if (ibpsw[0] == 1): |
|
for i in xrange(1, nrows - 1): |
|
W[i, 1] = 2.0 * W[i, 3 - 1] - W[i, 3] - pr * (W[i - 1, 3 - 1] - 2.0 * W[i, 3 - 1] + W[i + 1, 3 - 1]) |
|
|
|
for i in xrange(1, nrows - 1): |
|
comterm = 2.0 * (W[i, 3] - W[i, 1]) |
|
W[i, 0] = W[i, 4] - comterm + twompr * (W[i + 1, 3] - comterm + W[i - 1, 3] - W[i + 1, 1] - W[i - 1, 1]) |
|
|
|
else: |
|
for i in xrange(0, nrows): |
|
W[i, 0] = W[i, 3 - 1] |
|
|
|
|
|
# SET BOUNDARY CONDITIONS FOR RIGHT SIDE OF STUDY AREA |
|
if (ibpsw[1] == 1): |
|
for i in xrange(1, nrows - 1): |
|
W[i, ncols - 2] = 2.0 * W[i, ncols - 3] - W[i, ncols - 4] - pr * (W[i - 1, ncols - 3] - 2.0 * W[i, ncols - 3] + W[i + 1, ncols - 3]) |
|
|
|
for i in xrange(1, nrows - 1): |
|
comterm = 2.0 * (W[i, ncols - 4] - W[i, ncols - 2]) |
|
W[i, ncols - 1] = W[i, ncols - 5] - comterm + twompr * (W[i + 1, ncols - 4] - comterm + W[i - 1, ncols - 4] - W[i + 1, ncols - 2] - W[i - 1, ncols - 2]) |
|
|
|
else: |
|
for i in xrange(0, nrows): |
|
W[i, ncols - 1] = W[i, ncols - 3] |
|
|
|
|
|
# SET BOUNDARY CONDITIONS FOR BOTTOM OF STUDY AREA |
|
if (ibpsw[2] == 1): |
|
for j in xrange(1, ncols - 1): |
|
W[1, j] = 2.0 * W[3 - 1, j] - W[3, j] - pr * (W[3 - 1, j - 1] - 2.0 * W[3 - 1, j] + W[3 - 1, j + 1]) |
|
|
|
for j in xrange(1, ncols - 1): |
|
comterm = 2.0 * (W[4, j] - W[2, j]) |
|
W[0, j] = W[4, j] - comterm + twompr * (W[3, j + 1] - comterm + W[3, j - 1] - W[1, j + 1] + -W[1, j - 1]) |
|
|
|
else: |
|
for j in xrange(0, ncols): |
|
W[0, j] = W[3 - 1, j] |
|
|
|
|
|
# SET BOUNDARY CONDITIONS FOR TOP OF STUDY AREA |
|
if (ibpsw[3] == 1): |
|
for j in xrange(1, ncols - 1): |
|
W[nrows - 2, j] = 2.0 * W[nrows - 3, j] - W[nrows - 4, j] - pr * (W[nrows - 3, j - 1] - 2.0 * W[nrows - 3, j] + W[nrows - 3, j + 1]) |
|
|
|
for j in xrange(1, ncols - 1): |
|
comterm = 2.0 * (W[nrows - 3, j] - W[nrows - 1, j]) |
|
W[nrows - 1, j] = W[nrows - 5, j] - comterm + twompr * (W[nrows - 4, j + 1] - comterm + W[nrows - 4, j - 1] - W[nrows - 2, j + 1] - W[nrows - 2, j - 1]) |
|
|
|
else: |
|
for j in xrange(0, ncols): |
|
W[nrows - 1, j] = W[nrows - 3, j] |
|
|
|
|
|
# SET CORNER NODES TO AVERAGE OF 3 ADJACENT VALUES |
|
|
|
W[0, 0] = (W[0, 1] + W[1, 0] + W[1, 1]) / 3.0 |
|
W[0, ncols - 1] = (W[0, ncols - 2] + W[1, ncols - 1] + W[1, ncols - 2]) / 3.0 |
|
W[nrows - 1, 0] = (W[nrows - 1, 1] + W[nrows - 2, 0] + W[nrows - 2, 1]) / 3.0 |
|
W[nrows - 1, ncols - 1] = (W[nrows - 1, ncols - 2] + W[nrows - 2, ncols - 1] + W[nrows - 2, ncols - 2]) / 3.0 |
|
|
|
if not converged: |
|
print ' FAILED TO CONVERGE AFTER', its , 'ITERATION , MAX FRACTIONAL ERROR = ', change |
|
|
|
numpy.savetxt('DEF3D.DAT', W) |
|
return W |
|
|
|
def fortran_readline(F_obj): |
|
""" This will allow us to a read a line from the inumpy.t file in the same way fortran does |
|
USAGE example: |
|
PAR3DF=open('PAR3DF.DAT','rU') |
|
ym,pr = fortran_readline(PAR3DF) |
|
""" |
|
def s2n(s): |
|
if ('e' in s) or ('E' in s) or ('.' in s): |
|
n = float(s) |
|
else: |
|
n = int(s) |
|
return n |
|
|
|
line = F_obj.readline() |
|
splt = line.split(' ') |
|
content = splt[0] |
|
comment = splt[1:] |
|
nums = [s2n(x) for x in content.split(',')] |
|
if len(nums) == 1: nums = nums[0] |
|
return nums |
|
|
|
def fortran_readnrowncols(F_obj): |
|
""" This will allow us read ahead and give us the number of rows and cols. Yuck. |
|
""" |
|
ym, pr = fortran_readline(PAR3DF) |
|
ntv = fortran_readline(PAR3DF) |
|
for itv in xrange(0, ntv): |
|
ilo, ihi, jlo, jhi, thick = fortran_readline(PAR3DF) |
|
nl = fortran_readline(PAR3DF) |
|
for il in xrange(0, nl): |
|
ilo, ihi, jlo, jhi, sizeld = fortran_readline(PAR3DF) |
|
delx = fortran_readline(PAR3DF) |
|
nrows, ncols = fortran_readline(PAR3DF) |
|
return nrows, ncols |
|
|
|
if __name__ == "__main__": |
|
|
|
# This will allow us read ahead and give us the number of rows and cols |
|
PAR3DF = open('PAR3DF.DAT', 'rU') |
|
nrows, ncols = fortran_readnrowncols(PAR3DF) |
|
PAR3DF.close() |
|
|
|
# Load in settings file |
|
PAR3DF = open('PAR3DF.DAT', 'rU') |
|
ym, pr = fortran_readline(PAR3DF) |
|
print ym, pr |
|
|
|
# set up empty arrays |
|
W = numpy.zeros((nrows, ncols)) # DEFORMATION ARRAY |
|
D = numpy.zeros((nrows, ncols)) # FLEXURAL RIGIDITY ARRAY |
|
P = numpy.zeros((nrows, ncols)) # LOADING ARRAY |
|
GAMMA = numpy.zeros((nrows, ncols)) # BUOYANT RESTORING STRESS PER UNIT OF DISPLACEMENT |
|
IBPSW = numpy.zeros(4) # BROKEN PLATE SWITCH. FIRST ENTRY REFERS TO LEFT SIDE OF PLATE; SECOND TO RIGHT; THIRD TO BOTTOM; FOURTH TO TOP. "0" = CONTINUOUS; "1"=BROKEN. |
|
|
|
# calculate flexural RIGIDITY |
|
ntv = fortran_readline(PAR3DF) |
|
print ntv |
|
for itv in xrange(0, ntv): |
|
ilo, ihi, jlo, jhi, thick = fortran_readline(PAR3DF) |
|
print ilo, ihi, jlo, jhi, thick |
|
thick = thick * 1000.0 |
|
# CALCULATE FLEXURAL RIGIDITY AND STORE IN APPROPRIATE ENTRIES OF THE |
|
# ARRAY "D". |
|
fr = ym * thick ** 3 / (12.0 * (1 - pr ** 2)) |
|
for i in xrange(ilo - 1, ihi): |
|
for j in xrange(jlo - 1, jhi): |
|
D[i, j] = fr |
|
|
|
# READ IN LOADING PATTERN AND SAVE IN APPROPRIATE ENTRIES OF ARRAY "P". |
|
nl = fortran_readline(PAR3DF) |
|
print nl |
|
for il in xrange(0, nl): |
|
ilo, ihi, jlo, jhi, sizeld = fortran_readline(PAR3DF) |
|
print ilo, ihi, jlo, jhi, sizeld |
|
for i in xrange(ilo - 1, ihi): |
|
for j in xrange(jlo - 1, jhi): |
|
P[i, j] = sizeld |
|
|
|
# dexl= gridspacing |
|
delx = fortran_readline(PAR3DF) |
|
print delx |
|
delx = delx * 1000.0 |
|
del4 = delx ** 4 |
|
|
|
nrows, ncols = fortran_readline(PAR3DF) |
|
print nrows, ncols |
|
|
|
# READ IN BUOYANT RESTORING STRESS VALUES AND SAVE IN APPROPRIATE |
|
# ENTRIES OF ARRAY "GAMMA". |
|
ng = fortran_readline(PAR3DF) |
|
print ng |
|
for ig in xrange(0, ng): |
|
ilog, ihig, jlog, jhig, sizgam = fortran_readline(PAR3DF) |
|
print ilog, ihig, jlog, jhig, sizgam |
|
for i in xrange(ilog - 1, ihig): |
|
for j in xrange(jlog - 1, jhig): |
|
GAMMA[i, j] = sizgam |
|
|
|
|
|
# READ IN NO. OF ITERATIONS AND TOLERANCE FOR TESTING CONVERGENCE OF |
|
# GAUSS-SEIDEL ITERATION PROCEDURE. |
|
nits = fortran_readline(PAR3DF) |
|
print nits |
|
tol = fortran_readline(PAR3DF) |
|
print tol |
|
icl, ich, jcl, jch = fortran_readline(PAR3DF) |
|
print icl, ich, jcl, jch |
|
|
|
# "OMEGA" IS THE SIMULATANEOUS OVER-RELAXATION PARAMETER. GENERALLY |
|
# A VALUE OF 1.7 IS OK. |
|
omega = fortran_readline(PAR3DF) |
|
print omega |
|
ibpsw = fortran_readline(PAR3DF) |
|
print ibpsw |
|
fx, fy, fxy = fortran_readline(PAR3DF) |
|
print fx, fy, fxy |
|
|
|
# call the modelling program |
|
W = flex3db(D, P, GAMMA, delx, [0, 0, 0, 0], nits, tol) |
|
|
|
# basic visualisation |
|
try: |
|
import matplotlib.pyplot as plt |
|
except: |
|
print "No matplotlib installed so plotting will be skipped." |
|
else: |
|
plt.imshow(W, label = 'img') |
|
plt.colorbar() |
|
# plt.title('') |
|
# plt.xlabel('') |
|
# plt.ylabel('') |
|
plt.show() |
|
# plt.savefig('DEF3D.png') |
|
|
