Sunday, September 27, 2020

Cut and Fill Volume for given GL and FL profile, using Python Code

Hi,

Please don't get afraid by the length of the code. It is very simple to copy and paste it into your Python IDEs such as Pycharm or Visual Studio and run it. 


Once you run the code, it will ask you to input the section width and slope, and then the chainage, GL(ground level), FL(Floor level) etc. Output will be the cut or fill height(h), End Areas, and then the Cut and Fill Volume. The result will be in the tabular form. 

Code: 

import NumPy as np

from pandas import DataFrame

Chain = np.array([])
GL = np.array([])
FL = np.array([])
Cutting = np.array([])
Filling = np.array([])
Area = np.array([])
Vcut = np.array([0])
Vfill = np.array([0])
n = int(input('no. of sections, n = '))
b = float(input('width of level section b = '))
s = float(input('slope s = '))
cut = 0
fill = float()
A = float()

for i in range(n):
print(f'for section{i}')
ch = float(input('Chainage = '))
Chain = np.append(Chain, ch)
gl = float(input('GL = '))
GL = np.append(GL, gl)
fl = float(input('FL = '))
FL = np.append(FL, fl)
if gl>=fl:
fill = 0
Filling = np.append(Filling, fill)
cut = gl - fl
Cutting = np.append(Cutting, cut)
A = (b+s*cut)*cut
Area = np.append(Area, A)
else:
fill = fl - gl
Filling = np.append(Filling, fill)
cut = 0
Cutting = np.append(Cutting, cut)
A = (b + s * fill) * fill
Area = np.append(Area, A)
i+=1

print(f' Please note, When GL> FL, section is in cutting, h(cutting) = GL-FL,'
f'and if GL<FL it is in Filling, h(filling) = FL-GL'
f'Area of section, A = (b + sh)h', "\n")
df = DataFrame()
df['Chainage'] = Chain
df['GL'] = GL
df['FL'] = FL
df['Cutting'] = Cutting
df['Filling'] = Filling
df['Area'] = Area
print(df)

x = float()
D=float()
i=0
for i in range(n-1):

if Cutting[i]>0 and Filling[i+1]>0:
print(f'Both cutting and Filling in between Chainage {Chain[i]} to {Chain[i+1]}',)
D = Chain[i + 1] - Chain[i]
x = D*Cutting[i]/(Filling[i+1]+Cutting[i])
Vc = Area[i]*x/2
Vcut = np.append(Vcut,Vc)
Vf = Area[i+1]*(D-x)/2
Vfill = np.append(Vfill, Vf)
print(f'Cutting till distance x = {x}, and filling in {D-x}', "\n")

elif Filling[i]>0 and Cutting[i+1]>0:
print(f'Both cutting and Filling in between Chainage {Chain[i]} to {Chain[i + 1]}')
x = D * Filling[i] / (Filling[i] + Cutting[i+1])
print(f'Filling till distance x = {x}, and Cutting in {D - x}', "\n")
Vf = Area[i]*x/2
Vfill = np.append(Vfill, Vf)
Vc = Area[i+1]*(D-x)/2
Vcut = np.append(Vcut, Vc)

elif Cutting[i]>0 and Cutting[i+1]>0:
print(f'Only cutting in between Chainage {Chain[i]} to {Chain[i + 1]}')
Vc = D*(Area[i + 1] + Area[i]) / 2
Vcut = np.append(Vcut, Vc)
Vf = 0.0
Vfill = np.append(Vfill, Vf)
print(f' Cutting Volume = {D}*({Area[i]} + {Area[i+1]}) / 2 = {Vc}', "\n")

elif Filling[i]>0 and Filling[i+1]>0:
print(f'Only Filling in between chainage {Chain[i]} to {Chain[i + 1]}')
Vf = D * (Area[i + 1] + Area[i]) / 2
Vfill = np.append(Vfill, Vf)
Vc = 0.0
Vcut = np.append(Vcut, Vc)
print(f' Filling Volume = {D}*({Area[i]} + {Area[i + 1]}) / 2 = {Vc}')

i+=1

df['Filling Volume'] = Vfill
df['Cutting Volume'] = Vcut
print(df)

Friday, September 25, 2020

Python Code - Section Areas ( of Level Section, Two-Level Section, and Three-Level Sections ) and Cut and Fill Volume using Prismoidal Formula

 Hi,

When you have to find the Cut and Fill Volume of soil for a given profile, you have to find the section areas at a specified interval of chainage. End areas are calculated for a level trapezoid differently than for a two-level and three-level section. The following Python code contains all the formulas at the right places. 


Once you run the program, you have to put the inputs as per the standard formulas as and when asked. When asked for the unit, put a unit of length such as m, ft etc. Also, 'n' is for the number of sections, for example, you might need to calculate the areas at 5 different sections, on the given profile, then n = 5. 

The Volume Calculation at the end, for the cut and fill, is done using the prismoidal formula which works only when you have an odd number of sections. If you don't need the volume part, simply copy the code only above the comment "# Volume using Prismoidal Formula".

Code:

import math
import numpy as np
# Type of Sections
a = input('Section Areas given? y/n ')
A = float()
u = str()
l = str()
if a.lower() == "n":
l = input('Is there any level section? y/n ')

m = str()
if l.lower() == "n":
m = input('Is there any two-level section? y/n ')

w = str()
if l.lower() == "n":
w = input('Is there any three-level section? y/n ')
# Different approaches for level, two-level and three-level sections
Areas = np.array([])

if l.lower() == "y":
print("Level Trapezoidal Section,", "\n",
" area is given as A = (b+sh)h")
n = int(input('number of sections = '))
for i in range(n):
i+=1
print(f'for Secion{i}')
b = float(input('b = '))
s = float(input('s = '))
h = float(input('h = '))
u = input('unit= ')

A = (b+s*h)*h
Areas = np.append(Areas, A )
print(f' Area{i} A{i} = ({b}+{s}*{h})*{h} = {A} {u}^2', "\n")

if m.lower() == "y":
print("Two-Level Trapezoidal Section,", "\n",
" area is given as A = {(b/2s + h)(b1+b2)-b^2/2s}/2,", "\n"
"where, b1 = b/2 + n*s(h + b/(2*n))/(n-s)", "\n"
"b2 = b/2 + n*s(h - b/(2*n))/(n-s)", "\n")

n = int(input('number of sections = '))
for i in range(n):
i += 1
print(f'for Secion{i}')
b = float(input('b = '))
s = float(input('s = '))
h = float(input('h = '))
u = input('unit of length = ')
n = float(input('ground slope n = '))
b1 = b / 2 + n * s * (h + b / (2 * n)) / (n - s)
b2 = b / 2 + n * s * (h - b / (2 * n)) / (n - s)
A = ((b / 2 * s + h) * (b1 + b2) - b ** 2 / (2 * s)) / 2
Areas = np.append(Areas, A)

print(f''' putting all values for section{i}, b1 = {b1} {u}, b2= {b2} {u},
Area{i} A{i} = (({b} / 2 * {s} + {h}) * ({b1} + {b2}) - {b} ** 2 / (2 * {s})) / 2 '
f'= {A} {u}^2''', "\n")

if w.lower() == "y":

print("Three-Level Section,", "\n",
" area is given as A = h*(b1+b2)/2 + b*(h1+h2)/2", "\n"
"where, b1 = n1*s*(h + b/(2*s))/(n1-s)", "\n"
"b2 = n2*s*(h + b/(2*s))/(n2+s)"
"h1 = h + b1/n1"
" h2 = h - b2/n2", "\n")

n = int(input('number of sections = '))
for i in range(n):
i += 1
print(f'for secion{i}')
b = float(input('b = '))
s = float(input('s = '))
h = float(input('h = '))
u = input('unit= ')
n1 = float(input('ground slope1, n1 = '))
n2 = float(input('ground slope2, n2 = '))
b1 = n1 * s * (h + b / (2 * s)) / (n1 - s)
b2 = n2*s*(h + b/(2*s))/(n2+s)
h1 = h + b1/n1
h2 = h - b2/n2
A = h*(b1+b2)/2 + b*(h1+h2)/2
Areas = np.append(Areas, A)
print(f'putting all the values, b1= {b1}{u}, b2 = {b2}{u}, '
f'h1 = {h1}{u}, and h2 = {h2}{u}')
print(f' Area{i} A{i} = {h}*({b1}+{b2})/2 + {b}*({h1}+{h2})/2 '
f'= {A} {u}^2', "\n")

elif a.lower() == "y":
n = int(input('number of sections = '))
u = input('units of area = ')
for i in range(n):
i += 1
print(f'for Secion{i}')
A = input('Area = ')
Areas = np.append(Areas, A)

print(f'all the areas are {Areas} {u}^2')

# Volume using Prismoidal Formula

n = np.size(Areas) -1
print(f'number of areas, n = {n+1}')

V= float()
V2 = float()
V3 = float()
V1 = float()
if (n+1)%2 == 1:
print('Calculation of Cut and Fill using Prismoidal Formula:',
"\n", ' Volume (Cutting or Filling), V = D{A1 + An + 4(A2 + A4 + An-1) + 2(A3 + A5 + ..+An-2)}')
D = float(input('Distance between consecutive Sections, D = '))
V1 = D*(Areas[0] + Areas[n])/3
i = 0
while i < n/2:

V2 = V2+ 4*Areas[1+2*i]*D/3
i+=1
i=0
while i < (n/2-1):
V3 = V3 + 2 * Areas[2 + 2 * i] * D / 3
i+=1
V = V1 + V2 + V3
print(f'Putting all values, volume V = {V} {u}^3')

if (n+1)%2 == 0:
print('Prismoidal formular is not applicable for even number of sections')

Python Code for Drawing BMD and SFD (Bending Moment and Shear Force) for a SSB with Point Load

 Hi,

I have created this code in Python using Numpy and Matplotlib.pyplot libraries to draw the bending moment diagram(BMD) and Shear Force Diagram(SFD) for a Simply Supported Beam (SSB) with a point load at any location on that beam. 

The solution code when run will ask for the following inputs

Load, Load unit, Length of the beam, unit of the length, and location of the point load on the beam. 

And it will give the output for:

 1.Reactions,

2.Equation of the Shear and Moment at any location of the beam, 

3. Maximum Shear force and Maximum Bending Moment

4. Location of the maximum BM,

5. and Plot of the SFD and BMD.


Here is the code:

import numpy as np
import matplotlib.pyplot as plt

P = float(input('load = '))
u1 = input('load unit = ')
L = float(input('Length of the beam = '))
u2 = input('length unit = ')
a = float(input('Distance of Point load from left end = '))

b = L - a
R1 = P*b/L
R2 = P - R1
R1 = round(R1, 3)
R2 = round(R2, 3)
print(f'''
As per the static equilibrium, net moment sum at either end is zero,
hence Reaction R1 = P*b/L = {R1} {u1},
Also Net sum of vertical forces is zero,
hence R1+R2 = P, R2 = P - R1 = {R2} {u1}.
''')

l = np.linspace(0, L, 1000)
X = []
SF = []
M = []
maxBM= float()
for x in l:
if x <= a:
m = R1*x
sf = R1
elif x > a:
m = R1*x - P*(x-a)
sf = -R2
M.append(m)
X.append(x)
SF.append(sf)

print(f'''
Shear Force at x (x<{a}), Vx = R1 ={R1} {u1}
at x (x>{a}), SF = R1 - P = {R1} - {P} = -{R1-P} {u1}

Bending Moment at x (x<{a}), Mx = R1*x = {R1}*x
at x (x>={a}), Mx = R1*x - P*(x-{a})
= {R1}x - {P}(x-{a}) = -{R2}x + {P*a}
''')
max_SF = 0
for k in SF:
if max_SF < k:
max_SF = k


print(f'Maximum Shear Force Vmax = {max_SF} {u1}')


for k in M:
if maxBM < k:
maxBM = k
print(f'maximum BM, Mmax = {round(maxBM, 3)} {u1}{u2}')
Mx = float()

for x in l:
if x<a:
Mx = R1*x
if maxBM == Mx:
print(f'maximum BM at x = {round(x,3)} {u2}')
elif x>=a:
Mx = R1*x - P*(x- a)
if maxBM == Mx:
print(f'maximum BM at x = {round(x,3)} {u2}')


plt.plot(X, SF)
plt.plot([0, L], [0, 0])
plt.plot([0, 0], [0, R1], [L, L], [0, -R2])
plt.title("SFD")
plt.xlabel("Length in m")
plt.ylabel("Shear Force")
plt.show()


plt.plot(X, M)
plt.plot([0, L], [0, 0])
plt.title("BMD")
plt.xlabel("Length in m")
plt.ylabel("Bending Moment")
plt.show()

Monday, September 21, 2020

USCS and AASHTO- Soil Classification, using Software/Apps/Website building language Python

 Hi,

The following code is written in Python language, so to run the code you must have Python installed in your computer. Also, make sure to have 'math' library installed. 

When you run the code, it will ask you to enter the inputs related to the soil properties such as #200 sieve % finer, #4 finer etc. Put all the data given your problem, as and when asked. The result will show the USCS classification of the soil. 


def soilclass():
import math


r = input('Name of the classification required? (AASHTO/USCS/Both): ')
a = float(input('Material % finer than sieve #200 size(0.075mm): '))
b = float()
if a<50:
b = float(input('% finer than #4 (4.75mm) sieve: '))


LL = float()
PL = float()
PI = float()
PIA = float()

v = input('LL given? (y/n): ')
if v.lower() == "y":
LL = float(input('LL(%) = '))
pi = input('What is given? (PI or PL): ')
if pi.upper() == "PI":
PI= float(input('PI(%) = '))
PIA = 0.73 * (LL - 20)
if PI > PIA:
print(f'On A-Line, PI = 0.73(LL-20) = {PIA}%, '
f'Hence with PI of {PI}% soil lies above the A-line')
elif PI < PIA:
print(f'On A-Line, PI = 0.73(LL-20) = {PIA}%, '
f' Hence with PI of {PI}% soil lies below the A-line')
elif pi.upper() == "PL":
PL = float(input('PL(%)= '))
PI = LL - PL
print(f'Plasticity Index, PI = LL- PL = {PI}%')
PIA = 0.73 * (LL - 20)
if PI > PIA:
print(f'On A-Line, PI = 0.73(LL-20) = {PIA}%, '
f'Hence with PI of {PI}% soil lies above the A-line')
elif PI < PIA:
print(f'On A-Line, PI = 0.73(LL-20) = {PIA}%, '
f' Hence with PI of {PI}% soil lies below the A-line')

if r.upper() == "USCS" or r.upper() =="BOTH":
print('----------------')
print('USCS Classification: '
'----------')
D10 = float(input('D10(mm)= '))
D60 = float(input('D60(mm) = '))
D30 = float(input('D30(mm) = '))
Cu = D60 / D10
Cc = (D30) ** 2 / (D10 * D60)
print(f'''Cu = D60/D10 = {Cu}, Cc = (D30)^2/(D10*D60) = {Cc}''')

if a < 50:
print(f'As less than 50% of total soil, i.e. {a}%, is passing through #200 sieve, '
f'it is a Coarse grained soil')
if b < 50:
if a < 5:
if Cu >= 4 and 1 <= Cc <= 3:
print(f'As Cu>=4 and 1<=Cc<=3, '
f'soil class is GW (Well Graded Gravel)')
else:
print(f'As either Cu<4 or Cc<1 or Cc>3, '
f'soil class is GP (Poorly Graded Gravel)')
elif a > 12:
if PIA < PI and LL < 50:
print(f'With PI above A-line and LL<50, fines are clay with Low Compressibility(CH), '
f'so soil can be classified as GC - Clayey Gravel')
elif PIA < PI and LL > 50:
print(f'With PI above A-line and LL>50, fines are Clay with High Compressibility(CH), '
f'so soil can be classified as GC - Clayey Gravel')
elif PIA > PI and LL < 50:
print(f'With PI below A-line and LL<50,fines are Silt with Low Compressibility(ML), '
f'so soil can be classified as GM - Silty Gravel')
elif PIA > PI and LL > 50:
print(f'With PI below A-line and LL>50, fines are Silt with High Compressibility(MH), '
f'so soil can be classified as GC - Clayey Gravel')
exit()


elif b >= 50:
if a < 5:
if Cu >= 4 and 1 <= Cc <= 3:
print(f'As Cu>=4 and 1<=Cc<=3, '
f'soil class is SW (Well Graded Sand)')
else:
print(f'As either Cu<4 or Cc<1 or Cc>3, '
f'soil class is SP (Poorly Graded Sand)')

elif a > 12:
if PIA < PI and LL < 50:
print(f'With PI above A-line and LL<50, fines are clay with Low Compressibility(CH), '
f'so soil can be classified as SC (Clayey Sand)')
elif PIA < PI and LL > 50:
print(f'With PI above A-line and LL>50, fines are Clay with High Compressibility(CH), '
f'so soil can be classified as SC (Clayey Sand)')
elif PIA > PI and LL < 50:
print(f'With PI below A-line and LL<50, fines are Silt with Low Compressibility(ML), '
f'so soil can be classified as SM (Silty Sand)')
elif PIA > PI and LL > 50:
print(f'With PI below A-line and LL>50, fines are Silt with High Compressibility(MH), '
f'so soil can be classified as SC (Clayey Sand)')


elif a >= 50:
print(f'More than or equal to 50%, i.e. {a}%, is passing through #200 sieve')
if LL < 50:
if PIA < PI and PI > 7:
print(f'With PI above A-line and LL<50, '
f'fines are clay with Low Compressibility(CL)')
elif PIA > PI and PI < 4:
print(f'With PI below A-line, LL<50, and PI<4, '
f'fines are Silt with Low Compressibility(ML)')
elif LL > 50:

if PIA < PI:
print(f'With PI above A-line and LL>50, fines are Clay with High Compressibility(CH)')

elif PIA > PI:
print(f'With PI below A-line and LL>50, fines are Silt with High Compressibility(MH)', "\n")

if r.upper() == "AASHTO" or r.upper() == "BOTH":
print('---------------------- ')
print('AASHTO Soil Classification:')

# a is finer than #200, b is finer than #4
c = float(input('% finer than 2.00 mm (#10): '))
d = float(input('% finer than 0.425 mm (#40): '))

if c <= 50 and d <= 30 and a <= 15 and PI <= 6:
print(f'As % finer than 2.0mm (#200) <=50, %finer than 0.425 mm (#40)<=30'
f'% finer than 0.075mm (#200) <=15%, and PI<=6, '
f'hence the soil fits the criteria of soil class A-1-a (Stone Fragments,gravels and sand)')

elif c >= 51 and d <= 50 and a <= 25 and PI <= 6:
print(f'As %finer than 0.425 mm (#40)<=50'
f'% finer than 0.075mm (#200) <=25%, and PI<=6, '
f'hence the soil fits the criteria of soil class A-1-b (Stone Fragments,gravels and sand)')

elif c >= 51 and d >= 51 and a <= 10 and v.lower() == "n":
print(f'with %finer than 0.425 mm (#40) >=51, '
f'and % finer than 0.075mm (#200) <=10%, '
f'the soil is fine sand of class A-3')

elif c > 50 and d >=51 and a <= 35 and LL <= 40 and PI <= 10:
print(f'with LL<=40 and PI<=10, '
f'soil fits the criteria of "Silty or clayey Gravel" A-2-4')

elif c > 50 and d > 50 and a <= 35 and LL >= 41 and PI <= 10:
print(f'with LL>=41 and PI<=10, '
f'soil fits the criteria of "Silty or clayey Gravel" A-2-5')

elif c > 50 and d > 50 and a <= 35 and LL <= 40 and PI >= 11:
print(f'with LL<=40 and PI<=10, '
f'soil fits the criteria of "Silty or clayey Gravel" A-2-6')

elif c > 50 and d > 50 and a <= 35 and LL >= 41 and PI >= 11:
print(f'with LL>=41 and PI>=11, '
f'soil fits the criteria of "Silty or clayey Gravel" A-2-7')

elif c > 50 and d > 50 and a >= 36 and LL <= 40 and PI <= 10:
print(f'with LL<=40 and PI<=10, '
f'soil fits the criteria of "Silty Soils" A-4')

elif c > 50 and d > 50 and a >= 36 and LL >= 41 and PI <= 10:
print(f'with LL<=40 and PI<=10, '
f'soil fits the criteria of "Silty Soils" A-5')

elif c > 50 and d > 50 and a >= 36 and LL <= 40 and PI >= 11:
print(f'with LL<=40 and PI>=11, '
f'soil fits the criteria of "Silty Soils" A-6')

elif c > 50 and d > 50 and a >= 36 and LL >= 41 and PI >= 11:
if PI <= LL - 30:
print(f'with LL>=41 and PI>=11, and PI < LL-30 '
f'soil fits the criteria of "Silty Soils" A-7-5')
elif PI > (LL - 30):
print(f'with LL>=41 and PI>=11, and PI > LL-30 '
f'soil fits the criteria of "Silty Soils" A-7-5')


soilclass()

a = input('Exit? stress/n: ')
while a.lower() == "n":
soilclass()
a = input('Exit? stress/n: ')
else:
exit()

Sunday, September 20, 2020

Python Code for Mohr's Circle, - Stresse at a given angle theta, Principal stresses and their Planes, Maximum Shear Stress and its Plane

Hi,
You can use the Python to draw Mohr's Circle. In Structure Mechanics, taught to Civil Engineering and Mechanical
engineering students, one has to find the stresses at given oblique planes.
Principal stresses and maximum shear stresses are important to find, and also their 
corresponding planes. One can use the stress transformation equations or Mohr's circle method to 
find them. The following Python language code is prepared to be run on your IDE like Pycharm on your Computer.
One can also use online Jupyter's notebook to run the code. 
code: 
def mohcircle():
import numpy as np
import matplotlib.pyplot as plt
import math

σx = float(input('σx = '))
σy = float(input('σy = '))
τxy = float(input('τxy = '))
u = input('stress unit = ')
w = float(input("Angle( in degrees) of plane's axis from x axis(+ve counter clockwise), θ = "))
θ = math.radians(w)
R = np.sqrt(0.25 * (σx - σy) ** 2 + (τxy) ** 2)
σavg = (σx + σy) / 2
ψ = np.linspace(0, 2 * np.pi, 360)
x = σavg + R * np.cos(ψ)
y = R * (np.sin(ψ))
φ1 = math.degrees(0.5 * math.atan(2 * τxy / (σx - σy)))
φ2 = φ1 + 90
σθ1 = σavg + R * np.cos(2 * np.radians(φ1) + 2 * θ)
σθ2 = σavg + R * np.cos(2 * np.radians(φ1) + 2 * θ + np.pi)
τθ = R * np.sin(2 * np.radians(φ1) + 2 * θ)
print(f'''
Radius, R = √(0.25*(σx-σy)^2 + τxy^2)
= √(0.25*({σx}-{σy})^2 + {τxy}^2) ={R} {u}

Average Stress,(acts at the Center of Mohr's Circle)
= σavg = (σx + σy)/2 = ({σx} + {σy})/2 = {σavg} {u}

Principal Stresses
σ1 = σavg + R = {σavg} + {R} = {σavg + R} {u}
σ2 = σavg - R = {σavg} - {R} = {σavg - R} {u}
Angle φ1 it makes with x-axis,
φ1 = 0.5(atan(2*τxy/(σx - σy)) = 0.5 * atan(2*{τxy}/({σx} - {σy})) = {φ1} degrees
Angle σ2 makes with x-axis = φ2 = φ1 + 90 = {φ2} degrees

Maximum Shear Stress = τmax = R = {R} {u}
It occurs at, α = φ1 + 45 = {φ1 + 45} degrees

Stresses at a plane with axis at θ anticlockwise from x axis,
σθ1 = σavg + R* Cos(2φ1 + 2θ) = {σavg} + {R}* Cos({2 * φ1 + 2 * θ})
= {σθ1}, {u}
σθ2 = σavg + R* Cos(2φ1 + 2θ + pi) =
{σθ2} {u}
τθ = R*Sin(2*φ1 + 2*θ) = {R * np.sin(2 * np.radians(φ1) + 2 * θ)} {u}

''')

plt.plot(x, y)
plt.plot([σavg - R - 10, σavg + R + 10], [0, 0], linestyle='--', color='black')
plt.plot([σavg, σavg], [-R - 10, R + 10], linestyle='--', color='black')
plt.plot([σx, σy], [τxy, -τxy], [σx, σx], [0, τxy], [σy, σy], [0, -τxy], linestyle='-', color='brown')
plt.plot([σθ1, σθ2], [τθ, -τθ], [σθ1, σθ1], [0, τθ], [σθ2, σθ2], [0, -τθ], linestyle='--', color='red')
plt.xlabel('σ (MPa)')
plt.ylabel('τ (MPa)')
plt.title("Mohr's Circle")
plt.show()

mohcircle()
v = input('Exit? y/n ')
while v == "n":
mohcircle()
v = input('Exit? y/n ')
exit()

---- --------------
Thakns!

Friday, September 18, 2020

Solve Linear Equations in Three Variables - Using Python, Numpy and Sympy , With the Given Code

 Hi,

In the series of posts about Python for Civil Engineers, I have come with something from linear algebra.

While dealing with statics problems, such as finding unknown forces in space truss members, one has to deal with the system of linear equations in three variables. 


One can use the 'Sympy' and 'Numpy' libraries in Python can get a beautiful solution code for solving a system of linear equations in three variables. Computer programming is very much easy in Python, and you definitely need to try it. The code given below is separated into two parts: 1st part is the solution of the same problem using Sympy, and 2nd part gives solution code using Numpy library. 

The code is written specifically to solve the following three linear equations, but you can change the numbers accordingly.

x + 2y + 3z = 8

4x + 2y + 4z = 5

4x + 2y + 5z = 1

Simply change the coefficient for any other system. 


Code:

import numpy as np
import sympy as sp
from sympy import Eq, solve_linear_system, Matrix

# Using sympy

x,y,z = sp.symbols('x y z')
row1 = [1, 2, 3, 8]
row2 = [4, 2, 4, 5]
row3 = [4, 2, 5, 1]
system = Matrix((row1, row2, row3))
print(system)
print(solve_linear_system(system,x,y,z))

print('___________________________________ ')

# Using numpy
row1 = np.array([1,2,3])
row2 = np.array([4,2,4])
row3 = np.array([4,2,5])
constants = np.array([8,5,1])
system = np.array((row1, row2, row3))
print(system)
xyz= np.linalg.solve(system, constants)
print(xyz)


Thanks!

Thursday, September 17, 2020

Python Code to Find out the Dot Product of Vectors and Angle between them (Using end Point Co-ordinates too)

Hi, 

Sometimes you are given with the endpoints' co-ordinates of vectors, and sometimes you are directly given the vectors. For both conditions, the code given below will work. You have to input them as and when asked. End results will give the vectors, their magnitudes, their dot product and the angle between the two vectors. 

Python is a computer programming language that can be used to develip softwares, apps, and websites etc. The following python code provides a short solution to the Vector Algebra. 


Code: 


import numpy as np

import math

 a1 = input('For vector1, co-ordinates given? (y/n): ')

a2 = input('For vector2, co-ordinates given? (y/n): ')
u = input('unit: ')
if a1 == 'y':
x11= float(input('enter co-ordinates of starting point, x11= '))
y11 = float(input('y11= '))
z11 = float(input('z11 = '))
V11 = np.array([x11, y11, z11])

x21 = float( input('enter co-ordinates of end point, x21= '))
y21 = float(input('y21= '))
z21 = float(input('z21 = '))
V21 = np.array([x21, y21, z21])

print(f'Starting point, V11 = {V11}')
print(f'Ending point, V21 = {V21}')
v1 = V21 - V11
print(f'Vector1, V1 = V21 - V11 = {v1}')
v1m = np.sqrt((x21-x11) ** 2 + (y21-y11) ** 2 + (z21-z11) ** 2)
print(f" Magnitude of v1, v1m = {v1m} {u}")
else:
print('enter the elements in vector1 ')
x1 = float(input())
y1 = float(input())
z1 = float(input())
v1 = [x1,y1,z1]
print(f'v1 =({x1}i + {y1}j + {z1}k) {u}' '\n')
v1m = np.sqrt(x1**2 + y1**2 + z1**2)
print(f" Magnitude of v1, v1m = √(({x1})^2 + ({y1})^2 + ({z1})^2 = {v1m} {u}" '\n')

if a2 == 'y':
x12= float( input('enter co-ordinates of starting point, x21 = '))
y12 = float(input('y12= '))
z12 = float(input('z12 = '))
V12 = np.array([x12, y12, z12])

x22 = float( input('enter co-ordinates of end point, x22 = '))
y22 = float(input('y22= '))
z22 = float(input('z22 = '))
V22 = np.array([x22, y22, z22])
print(f'Starting point, V12 = {V12}')
print(f'Ending point, V22 = {V22}')
v2 = V22 - V12
print(f'Vector2, V2 = V22 - V12 = {v2}')
v2m = np.sqrt((x22 - x12) ** 2 + (y22 - y12) ** 2 + (z22 - z12) ** 2)
print(f" Magnitude of v2, v2m = {v2m} {u}")

else:
print('enter the elements in the 2nd vector: ')
x2 = float(input())
y2 = float(input())
z2 = float(input())
v2 = [x2,y2,z2]
print(f'v2 = ({x2}i + {y2}j + {z2}k) {u}' '\n')
v2m = np.sqrt(x2**2 + y2**2 + z2**2)
print(f" Magnitude of v2, v2m = √(({x2})^2 + ({y2})^2 + ({z2})^2 = {v2m} {u}" '\n')

d = float(np.dot(v1,v2))
e = d/(v1m*v2m)
a = math.acos(e)
k = math.degrees(a)

print(f'''Taking Dot product, v1.v2 = v1m * v2m * cosθ

=> cosθ = v1.v2/(v1m*v2m) = {d}/({v1m}*{v2m}) = {e}

hence, Angle between the two vectors, θ = acos({e}) = {k}
''')

Cut and Fill Volume for given GL and FL profile, using Python Code

Hi, Please don't get afraid by the length of the code. It is very simple to copy and paste it into your Python IDEs such as Pycharm or V...