# Using SciPy and Other Third-party Python Packages

This example demonstrates using an optimization routine from `scipy`

to solve an inverse elastic problem.

A constant force, f, is applied to the tip of a cantilever beam. From a given beam geometry and displacement we use PFC to find the modulus of the beam. The beam is modeled with parallel bonded PFC particles.

The deflection is known analytically,

where f is the applied force, E is the Youngs Modulus, I is the moment of inertia and x is the length along the beam and L is the total length of the beam.

This is a simplification of a PFC verification problem see Cantilever Beam for more information.

```
import matplotlib
matplotlib.use("QT5Agg")
matplotlib.rcParams['backend.qt5']='PySide2'
import pylab as plt
import numpy as np
from scipy.constants import inch
import itasca as it
it.command("python-reset-state false")
from itasca import ballarray as ba
r = 8.0 # beam radius [inches]
l = 200.0 # beam length [inches]
x = np.linspace(0, l, 100)
E = 1e4 # [ksi]
I = 0.25*np.pi*r**4 # moment of inertia.
f = 1.0 # [kip] force applied in the positive y direction
```

First we import the required libraries and setup some parameters for the beam.

```
def analytical_deflection(x):
"Return the analytical solution for y deflection at the x point(s)"
return f * x**2 / 2.0 / E / I * ( l - x / 3.0)
def pfc_deflection(E_trial):
"""For a given parallel bond stiffness, return a tuple containing the
ball x positions and the y displacements.
calling this function cycles the PFC engine
"""
n = 11 # number of pfc particles
pb_kn = E_trial/l*(n-1)
print ("calculating deflection with pb_kn = {} ... ".format(pb_kn),)
it.command("""
model new
model large-strain on
model random
model domain ext -50 250 -50 50
contact cmat default model linearpbond
""")
for i in range(n):
it.ball.create(10.0, (i*20.0, 0.0, 0.0))
tip_ball = it.ball.find(n)
tip_ball.set_force_app_y(f)
it.command("""
ball attribute dens=1000 damp 0.7
ball property 'kn'=0 'ks'=0
model clean
contact method bond
contact property pb_ten 1e5 pb_coh 1e5 pb_kn {kn} pb_ks 188 pb_rmul 0.8
ball fix velocity-x velocity-y velocity-z ...
spin-x spin-y spin-z range id=1
model solve ratio-maximum 1e-4
""".format(kn=pb_kn))
pfc_x = ba.pos()[:,0]
pfc_dy = ba.disp()[:,1]
print ("done")
return pfc_x, pfc_dy
```

Two functions are defined one to give the analytical solution and one to give the PFC solution.

We can test everything is working by calling the pfc_deflection function with the known correct modulus. The pfc and analytical solution should match,

```
pfc_x, pfc_dy = pfc_deflection(E)
plt.plot(x, analytical_deflection(x))
plt.plot(pfc_x, pfc_dy, "o")
plt.xlabel("Length along beam [inches]")
plt.ylabel("Beam vertical displacement [inches]")
plt.draw()
plt.savefig("p3d-opt_comp.png")
```

- output:
calculating deflection with pb_kn = 500.0 ... done

The PFC beam deflection is shown as points. The PFC solution matches well against the analytical solution which is shown as a solid blue line.

Next we do the actual optimization calculation. The function
`fminbound`

is imported from the `scipy.optimize`

module. This
function takes three arguments, a callable function, an upper bound
and a lower bound. We define the function `residual`

to return the
sum of the square of the difference between the PFC and analytical
solutions. We bound the optimization between 1e2 and 5.5e4.

- For more about
`fminbound`

see: http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html

```
from scipy.optimize import fminbound
def residual(pb_kn):
"""Return the sum of squares of the residual. This is a scalar measure
of the difference between the PFC solution and the analytical solution.
"""
pfc_x, pfc_dy = pfc_deflection(pb_kn)
a_dy = analytical_deflection(pfc_x)
plt.close('all')
plt.plot(x, analytical_deflection(x))
plt.plot(pfc_x, pfc_dy, "o")
plt.xlabel("Length along beam [inches]")
plt.ylabel("Beam vertical displacement [inches]")
plt.draw()
return ((a_dy - pfc_dy)**2).sum()
E_opt = fminbound(residual, 1e2, 5.5e4)
print ("""
E specified: {}
E calculated: {}""".format(E, E_opt))
plt.close('all')
```

- output:
calculating deflection with pb_kn = 1053.4967008815386 ... done calculating deflection with pb_kn = 1701.5032991184614 ... done calculating deflection with pb_kn = 653.0065982369227 ... done calculating deflection with pb_kn = 405.490102644616 ... done calculating deflection with pb_kn = 530.7543627593386 ... done calculating deflection with pb_kn = 529.6311596321268 ... done calculating deflection with pb_kn = 482.21349526223526 ... done calculating deflection with pb_kn = 501.8190835557665 ... done calculating deflection with pb_kn = 501.70523214550855 ... done calculating deflection with pb_kn = 499.7989842280888 ... done calculating deflection with pb_kn = 493.08192515191905 ... done calculating deflection with pb_kn = 500.68789623953967 ... done calculating deflection with pb_kn = 500.5925746896699 ... done calculating deflection with pb_kn = 500.66760871604833 ... done calculating deflection with pb_kn = 500.6732746302266 ... done calculating deflection with pb_kn = 500.6776216232916 ... done calculating deflection with pb_kn = 500.6759566385157 ... done calculating deflection with pb_kn = 500.681546177477 ... done calculating deflection with pb_kn = 500.67693683977063 ... done calculating deflection with pb_kn = 500.67762921620755 ... done calculating deflection with pb_kn = 500.6791253622799 ... done calculating deflection with pb_kn = 500.678200693155 ... done calculating deflection with pb_kn = 500.6778475009777 ... done calculating deflection with pb_kn = 500.6777125935705 ... done calculating deflection with pb_kn = 500.6776610635263 ... done calculating deflection with pb_kn = 500.6776413808009 ... done E specified: 10000.0 E calculated: 10013.55258432415

The calculated value for E and the known value are close.

## Example File

The source code for the examples here is using_scipy.py.

Was this helpful? ... | 3DEC © 2019, Itasca | Updated: Feb 25, 2024 |