The source code for this example is using_scipy.py

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,

\[\Delta y = \frac {f x^{2}}{2EI} \left(L - \frac{x}{3}\right)\]

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 random
    model domain ext -50 250 -50 50
    model 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.

pfcpythonmodule/test3d/examples/test_ucs/simple/p3d-opt_comp.*

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.

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.9960911725644 ...
done
calculating deflection with pb_kn = 529.6933321485365 ...
done
calculating deflection with pb_kn = 482.2519199905426 ...
done
calculating deflection with pb_kn = 501.7128466491913 ...
done
calculating deflection with pb_kn = 501.7500039455624 ...
done
calculating deflection with pb_kn = 500.3288224273459 ...
done
calculating deflection with pb_kn = 493.42406010780286 ...
done
calculating deflection with pb_kn = 500.8378391908291 ...
done
calculating deflection with pb_kn = 500.67206251922556 ...
done
calculating deflection with pb_kn = 500.64970317429237 ...
done
calculating deflection with pb_kn = 500.73522749612437 ...
done
calculating deflection with pb_kn = 500.69618939350227 ...
done
calculating deflection with pb_kn = 500.68570703670133 ...
done
calculating deflection with pb_kn = 500.71110062184766 ...
done
calculating deflection with pb_kn = 500.6924573898482 ...
done
calculating deflection with pb_kn = 500.6898789843822 ...
done
calculating deflection with pb_kn = 500.6932356808154 ...
done
calculating deflection with pb_kn = 500.6924497967122 ...
done
calculating deflection with pb_kn = 500.6927546705446 ...
done
calculating deflection with pb_kn = 500.69257094097003 ...
done
calculating deflection with pb_kn = 500.6925007625173 ...
done
calculating deflection with pb_kn = 500.6924739567336 ...
done
calculating deflection with pb_kn = 500.6924649829842 ...
done


E specified: 10000.0
E calculated: 10013.849147796964

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