<jupyter_start><jupyter_text>Finite Element Analysis (1-D Elements)   <NAME>  @iampramodyadav ---![Imgur](https://i.imgur.com/nTr2w2q.jpg "FEA integration Flow Chart")  Natural Coordinate: Shape Funtion and its derivative<NAME> (@iampramodyadav) --- - This programme calculates the values of Shape functions and their derivatives.- plot the Shape funtions  Function for Shape function- [Nptel Reference](https://nptel.ac.in/content/storage2/courses/112104116/lecture-11/11_10.htm)- [SymPy Documentation](https://docs.sympy.org/latest/tutorial/basic_operations.html)- **Shape Functions in Natural Coordinate**  - ![](https://i.imgur.com/1iBlgQQ.png)  -Here,Approximation of Order p,   - physical coordinate ($x$) and natural coordinate $\xi$ (here, z)<jupyter_code>from sympy import *
import numpy as np
z=Symbol('z')
def Shape(p):
  '''
  =====================================================
  Shape(p)
  p:  order (p) of approximation
  z:  value of natural coordinate
  This function return Shape funtions values at given p
  =====================================================
  '''
  z=Symbol('z')
  n=[]

  for i in range(0, p+1):
    point=-1
    point=point+2*i/p
    n.append(point)

  shape=[1]*(p+1)

  for i in range(0,p+1):
    for j in range(0,p+1):
      if i!=j:
        shape[i]=shape[i]*((z-n[j])/(n[i]-n[j]))

  return shape
print(Shape.__doc__)
Shape(2)<jupyter_output><empty_output><jupyter_text>---  Function for derivative of Shape function- [SymPy Documentation](https://docs.sympy.org/latest/tutorial/calculus.html)- [Reference](https://stackoverflow.com/questions/47998930/differentiation-of-a-list-of-functions-with-respect-to-a-list-of-variables)<jupyter_code>def dShape(p):
  '''
  =====================================================
  dShape(p)
  p:  order (p) of approximation

  z:  natural coordinate variable
  This function return derivative of Shape funtions values at given p
  =====================================================
  '''  
  z=Symbol('z')
  f=Shape(p)
  z=symbols('z')
  return [simplify(diff(f1)) for f1 in f]
print(dShape.__doc__)<jupyter_output>=====================================================
  dShape(p)
  p:  order (p) of approximation

  z:  natural coordinate variable
  This function return derivative of Shape funtions values at given p
  =====================================================<jupyter_text>---  Examples<jupyter_code>df=dShape(3)
print('Derivative of Shape function:\n',df)
f=Shape(3)
print('Shape funtion:\n',f ,'\n')

print('Simplified Shape funtion:\n',[simplify(f1) for f1 in f])
simplify(f[0]) #first term
df[0] #deruvative of first term<jupyter_output><empty_output><jupyter_text>---  Visualization: Plotting Shape functions<jupyter_code>import seaborn as sns
def plotShape(p):
  '''
  =====================================================
  plotShape(p)

  p: order (p) of approximation
  z:value of natural coordinate
  This function return plot of Shape funtions (Natural Coordinate)
  =====================================================
  '''
  sns.set()
  sns.set_style("whitegrid", {'grid.linestyle': '--'})

  z=Symbol('z')
  f=Shape(p)
  p1=plot(0,(z,-1,1),show=False)
  for i in f:
    p2=plot(i,(z,-1,1),show=False)
    p1.append(p2[0])

  return p1.show()
print(plotShape.__doc__)
plotShape(1)
plotShape(2)
plotShape(3)
plotShape(4)
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 14, 8

plotShape(5)<jupyter_output><empty_output><jupyter_text>---  Gauss-Legendre Quadrature implementation in PythonPramod <NAME> (@iampramodyadav) ---   Important Link and Formula Used* [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature)* [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomialsRecurrence_relations)* [Newton–Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method)* **Gaussian quadrature Rule**$${\displaystyle \int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i})}$$* **Bonnet’s recursion formula for Legendre Polynomial**    * $$\boxed{\displaystyle (n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x)}$$    * $${\displaystyle P_{0}(x)=1\,}$$    * $${\displaystyle\quad P_{1}(x)=x}$$* Recursive equation for derivative of Legendre Polynomial:  * $$\boxed{\displaystyle {\frac {x^{2}-1}{n}}{\frac {d}{dx}}P_{n}(x)=xP_{n}(x)-P_{n-1}(x)\,}$$* The $i-th$ Gauss node, $x_i$, is the $i-th$ root of $P_n$ and the weights are given by the formula (Abramowitz & Stegun 1972)$$\boxed{\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.}$$* Newton–Raphson method$${\displaystyle x_{1}=x_{0}-{\frac {f(x_{0})}{f'(x_{0})}}}$$* The first guess ${\displaystyle x_{0}}$ for the ${\displaystyle i}$-th root of a ${\displaystyle n}$-order polynomial ${\displaystyle P_{n}}$ can be given by$${\displaystyle x_{0}=\cos \left({\frac {\pi (i-1/4)}{n+1/2}}\right)}$$  Legendre Polynomial (symbolic)[SymPy Documentation](https://docs.sympy.org/latest/tutorial/basic_operations.html)**Bonnet’s recursion formula for Legendre Polynomial**  *  $$\boxed{\displaystyle (n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x)}$$      * $${\displaystyle P_{0}(x)=1\,}$$  * $${\displaystyle\quad P_{1}(x)=x}$$<jupyter_code>from sympy import *
import numpy as np
x=symbols('x')
def Legendre(n):
  """
  =====================================================
  n: Order of polynomial
  x: Variable
  This function print Legendre polynomial of order n (Symbolic)
  =====================================================
  """

  x=symbols('x')
  if (n==0):
    return x*0+1.0
  elif (n==1):
    return x
  else:
    return ((2.0*n-1.0)*x*Legendre(n-1)-(n-1)*Legendre(n-2))/n
print(Legendre.__doc__)<jupyter_output>=====================================================
  n: Order of polynomial
  x: Variable
  This function print Legendre polynomial of order n (Symbolic)
  =====================================================<jupyter_text>**Example:**<jupyter_code>Legendre(2)
Legendre(5)
simplify(Legendre(5))<jupyter_output><empty_output><jupyter_text>Derivative of the Legendre polynomials[Reference](https://en.wikipedia.org/wiki/Legendre_polynomialsRecurrence_relations)$$\boxed{\displaystyle {\frac {x^{2}-1}{n}}{\frac {d}{dx}}P_{n}(x)=xP_{n}(x)-P_{n-1}(x)\,}$$<jupyter_code>def dLegendre(n):
  """
  =====================================================
  n: Order of polynomial
  x: Variable
  This function print Derivative of Legendre polynomial of order n  (Symbolic)
  =====================================================
  """
  x=symbols('x')
  if (n==0):
    return x*0
  elif (n==1):
    return x*0+1.0
  else:
    return (n/(x**2-1.0))*(x*Legendre(n)-Legendre(n-1))
print(dLegendre.__doc__)<jupyter_output>=====================================================
  n: Order of polynomial
  x: Variable
  This function print Derivative of Legendre polynomial of order n  (Symbolic)
  =====================================================<jupyter_text>**Example:**<jupyter_code>simplify(dLegendre(2))
dLegendre(5)
simplify(dLegendre(5))<jupyter_output><empty_output><jupyter_text>Visualization: Plot Legendre Polynomial<jupyter_code>f=Legendre(2)
from matplotlib import style
import matplotlib.pyplot as plt

#plt.rcParams['figure.figsize'] = 14, 8
style.use('ggplot')
plot(f,(x,-1,1),legend=True)

f1=Legendre(1)
f2=simplify(Legendre(3))
pl=plot(f1,f2,(x,-1,1),legend=True,show=False)
pl[0].line_color = 'red'
pl[1].line_color = 'green'
pl.show()
plot(Legendre(1),Legendre(2),Legendre(3),Legendre(4),Legendre(5),(x,-1,1))
import seaborn as sns
sns.set()
plt.rcParams['figure.figsize'] = 14, 8
sns.set_style("whitegrid", {'grid.linestyle': '--'})
p=plot(Legendre(1),Legendre(2),Legendre(3),Legendre(4),Legendre(5),(x,-1,1),title='Legendre PolyNomial',show=False);
p[0].line_color = 'red'
p[1].line_color = 'green'
p[2].line_color = 'yellow'
p[3].line_color = ['b','r']

p.show()<jupyter_output><empty_output><jupyter_text>Legedre Roots* The first guess ${\displaystyle x_{0}}$ for the ${\displaystyle i}$-th root of a ${\displaystyle n}$-order polynomial ${\displaystyle P_{n}}$ can be given by$${\displaystyle x_{0}=\cos \left({\frac {\pi (i-1/4)}{n+1/2}}\right)}$$* Newton–Raphson method$${\displaystyle x_{1}=x_{0}-{\frac {f(x_{0})}{f'(x_{0})}}}$$<jupyter_code>def LegendreRoots(polyorder,tolerance=1e-20):
  """
  =====================================================
  LegendreRoots(polyorder,tolerance=1e-20)

  polyorder: Order of polynomial(>2)
  tolerance: tolerance of error (default: 1e-20)
  Output: [roots, err]
  This function return Roots of Legendre polynomial
  =====================================================
  """
  if polyorder<2:
    err=1
  
  else:
    roots=[]
		                                                  # The polynomials are alternately even and odd functions. 
                                                      # So we evaluate only half the number of roots. 
    for i in range(1,int((polyorder)/2) +1):
      x1=cos(pi*(i-0.25)/(polyorder+0.5))
      error=10*tolerance
      iters=0

      while (error>tolerance) and (iters<1000):
        dx=-Legendre(polyorder)/dLegendre(polyorder)
        dx=N(dx.subs(x,x1))
        x1=N(x1+dx)
        iters=iters+1
        error=abs(dx)

      #print(roots)
      roots.append(x1)
      #print(roots)
		                                                  # Use symmetry to get the other roots
    roots=np.array(roots)
    if polyorder%2==0:
      roots=np.concatenate( (-1.0*roots, roots[::-1]) )

    else:
      roots=np.concatenate( (-1.0*roots, [0.0], roots[::-1]) )

    err=0
                                                    # successfully determined roots
  return [roots, err]
print(LegendreRoots.__doc__)<jupyter_output>=====================================================
  LegendreRoots(polyorder,tolerance=1e-20)

  polyorder: Order of polynomial(>2)
  tolerance: tolerance of error (default: 1e-20)
  Output: [roots, err]
  This function return Roots of Legendre polynomial
  =====================================================<jupyter_text>**Examples**<jupyter_code>[Xis, err]=LegendreRoots(8)
print(Xis)
[Xis, err]=LegendreRoots(2)
print(Xis)
[Xis, err]=LegendreRoots(4)
print(Xis)
[Xis, err]=LegendreRoots(5)
print(Xis)<jupyter_output>[-0.906179845938664 -0.538469310105683 0.0 0.538469310105683
 0.906179845938664]<jupyter_text>Weight Coefficient* the $i-th$ Gauss node, $x_i$, is the $i-th$ root of $P_n$ and the weights are given by the formula (Abramowitz & Stegun 1972)$${\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.}$$<jupyter_code>def GaussLegendreWeights(polyorder):
  """
  =====================================================
  GaussLegendreWeights(polyorder)

  polyorder: Order of polynomial
  Output: [W, xis, err]
  This function return Weights(W) and Roots(xis) of Legendre polynomial
  =====================================================
  """
  W=[]
  [xis,err]=LegendreRoots(polyorder)
  xis=list(xis)
  if err==0:
    for x1 in xis:
      w=2.0/( (1.0-x**2)*(dLegendre(polyorder)**2) )
      w=w.subs(x,x1)
      W.append(w)
      err=0
  else:
    err=1 # could not determine roots - so no weights
  return [W, xis, err]
print(GaussLegendreWeights.__doc__)<jupyter_output>=====================================================
  GaussLegendreWeights(polyorder)

  polyorder: Order of polynomial
  Output: [W, xis, err]
  This function return Weights(W) and Roots(xis) of Legendre polynomial
  =====================================================<jupyter_text>**Example**<jupyter_code>[W, xis, err]=GaussLegendreWeights(2)
print(W)
[W, xis, err]=GaussLegendreWeights(5)
print(W)<jupyter_output>[0.236926885056189, 0.478628670499366, 0.568888888888889, 0.478628670499366, 0.236926885056189]<jupyter_text>Function Plot Legendre Polynomial<jupyter_code>def PlotLegendre(polyorder):
  """
  =====================================================
  PlotLegendre(polyorder)

  This function plot legendre plynomial
  =====================================================
  """
  sns.set()
  sns.set_style("whitegrid", {'grid.linestyle': '--'})
  plot(Legendre(polyorder),(x,-1,1))

PlotLegendre(4)
PlotLegendre(2)<jupyter_output><empty_output><jupyter_text>Testing Final Result<jupyter_code>for Polyorder in range(2,7):
  plt.rcParams['figure.figsize'] = 6.4, 4.8
  plot(Legendre(Polyorder),(x,-1,1))

  [W, xis, err]=GaussLegendreWeights(Polyorder)
  print("order:", Polyorder)
  print("Weights:", W)
  print ("Roots    : ",xis)

  print('='*110)
  print('\n\n')
  print('='*110)<jupyter_output><empty_output><jupyter_text>---  Values of Shape Function at Gauss Points (Legendre root)<NAME> (iampramodyadav) Flow chart for Gauss-Legendre Quadrature use in FEA![Imgur](https://i.imgur.com/nTr2w2q.jpg "FEA integration Flow Chart") ---<jupyter_code>def ShapeLegRoot(p):
  """
  =====================================================
  ShapeLegRoot(p)

  p : order of approximation (integer)
  
  Output: [n, W, xis, ShapeValueXi, dShapeValueXi]
  Here,
  n-    Order of gauss legendre polynomial
  W-    Legendre Weights
  xis-  Gauss point coordinates (real)
  ShapeValueXi-   Value of Shape function at  at all the  Xi (real)(Legendre Roots)
  ShapeValueXi-   Value of derivative of Shape function at Gauss points Xi(Legendre Roots) 
  =====================================================
  """

  p=int(p)

  pbar=max(2*p, p+2) 
  if pbar%2!=0:                                   #if odd
    n=(pbar+1)/2
  else:
    n=(pbar+2)/2

  n=int(n)

  #print("Order of gauss legendre polynomial:",n)

  [W, xis, err]=GaussLegendreWeights(n)           #order: Polyorder(n), Weights:  W, Roots: xis

  ShapeValueXi=np.zeros((n, p+1))
  dShapeValueXi=np.zeros((n, p+1))

  #print('ndofel(p+1):',p+1)
  #print('Shape of Shape value matrix:',ShapeValueXi.Shape)

  s=Shape(p)
  ds=dShape(p)

  for i in range(0,n):
    for j in range(0,p+1):

      ShapeValueXi[i][j]=s[j].subs(z,xis[i])          #w.subs(x,x1)
      dShapeValueXi[i][j]=ds[j].subs(z,xis[i])

  return [n, W, xis, ShapeValueXi, dShapeValueXi]
print(ShapeLegRoot.__doc__)<jupyter_output>=====================================================
  ShapeLegRoot(p)

  p : order of approximation (integer)
  
  Output: [n, W, xis, ShapeValueXi, dShapeValueXi]
  Here,
  n-    Order of gauss legendre polynomial
  W-    Legendre Weights
  xis-  Gauss point coordinates (real)
  ShapeValueXi-   Value of Shape function at  at all the  Xi (real)(Legendre Roots)
  ShapeValueXi-   Value of derivative of Shape function at Gauss points Xi(Legendre Roots) 
  =====================================================<jupyter_text>--- **Example**<jupyter_code>[n, W,xis,ShapeValueXi,dShapeValueXi]=ShapeLegRoot(1)
print('Shape Function values at Legendre roots:\n',ShapeValueXi)
print('Derivative of Shape Function values at Legendre roots\n:',dShapeValueXi)
[n, W,xis,ShapeValueXi,dShapeValueXi]=ShapeLegRoot(5)
print('Shape Function values at Legendre roots:\n',ShapeValueXi)
print('Derivative of Shape Function values at Legendre roots\n:',dShapeValueXi)<jupyter_output>Shape Function values at Legendre roots:
 [[ 0.66464722  0.6750085  -0.61277567  0.39633694 -0.14644321  0.02322621]
 [ 0.04145294  1.14719867 -0.30450086  0.16307144 -0.05567618  0.008454  ]
 [-0.00517461  0.05451099  1.0201782  -0.08982382  0.02349007 -0.00318084]
 [-0.00318084  0.02349007 -0.08982382  1.0201782   0.05451099 -0.00517461]
 [ 0.008454   -0.05567618  0.16307144 -0.30450086  1.14719867  0.04145294]
 [ 0.02322621 -0.14644321  0.39633694 -0.61277567  0.6750085   0.66464722]]
Derivative of Shape Function values at Legendre roots
: [[-4.2710757   7.68824036 -5.97290352  3.67209091 -1.32305421  0.20670216]
 [-0.87306492 -2.03347341  4.8542633  -2.76385889  0.96414612 -0.14801221]
 [ 0.14181715 -1.57319431 -0.20311928  2.13898231 -0.58491528  0.08042941]
 [-0.08042941  0.58491528 -2.13898231  0.20311928  1.57319431 -0.14181715]
 [ 0.14801221 -0.96414612  2.76385889 -4.8542633   2.03347341  0.87306492]
 [-0.20670216  1.32305421 -3.67209091  5.97290352 -7.68824036  4.2710757 ]]<jupyter_text>Element Calculation * Here linear mapping function used between the physical coordinate (x) and the natural coordinate $\xi$ (here, z) :  * ![](https://i.imgur.com/2yMRAeP.png)  * Hence, $x_1^k$ (xL) and $x_{p+1}^k$ (xR) are the coordinates of the end nodes of the element $k$ and $p$ denotes the order of approximation.* Expressions for Element Stiffness Matrix and Force Vector  * $$\boxed{\displaystyle k_{i j}^{k}=\int_{-1}^{1} \hat{E}(\xi) \hat{A}(\xi) \frac{d \hat{N}_{i}(\xi)}{d \xi} \frac{d \hat{N}_{j}(\xi)}{d \xi} \frac{2}{h_{k}} d \xi \quad \text { for } i, j=1,2, \ldots, p+1}$$  * $$\boxed{\displaystyle f_{i}^{k}=\int_{-1}^{1} \hat{f}(\xi) \hat{N}_{i}(\xi) \frac{h_{k}}{2} d \xi \quad i=1,2, \ldots, p+1}$$    Function for Local $K_{ij}$, & $f_i$ if 1-D Bar Elements<jupyter_code>def ElemKfi(k,p,xL,xR,E,a0,a1=0,a2=0,f0=0,f1=0,f2=0):
  """
  =====================================================
  ElemKfi(k,p,xL,xR,E,a0,a1=0,a2=0,f0=0,f1=0,f2=0)

  k :   element number (integer)
  p :   order of approximation (integer)
  xL :  x -coordinate of the left end of the element (real)
  xR :  x -coordinate of the right end of the element (real)
  E :   Young's modules ( E ) of the bar material which is taken as constant (real)
  a0,a1,a2 : coefficients in (A=a0 + a1*x + a2*(x**2)) for the area of cross-section ( A ) of the bar (real)
  f0,f1,f2 : coefficients in (f = f0 + f1*x + f2*x*x) for the distributed force ( f ) acting on the bar (real)

  Output: [Ke, fe]
  Ke : element stiffness matrix (real)
  fe : element force vector (real)
  =====================================================
  """

  Ke=np.zeros((p+1,p+1))  # size of element matrix
  h=xR - xL               # h : element length (real)

  [n, W,xis,ShapeValueXi,dShapeValueXi]=ShapeLegRoot(p)
  fe=np.zeros(p+1)

  for m in range(0,n):
    xi=xis[m]
    x = xL*(1 - z)/2 + xR*(1 + z)/2

    EA = E*(a0 + a1*x+a2*(x**2))
    EA = EA.subs(z,xi)

    f = f0 + f1*x + f2*x*x
    f=f.subs(z,xi)

    for i in range(0, p+1):
      fe[i] = fe[i] + f*ShapeValueXi[m][i]*(h/2)*W[m]

      for j in range(0,p+1):
        Ke[i][j] = Ke[i][j]+EA*dShapeValueXi[m][i]*dShapeValueXi[m][j]*(2/h)*W[m]

  return [Ke, fe]
print(ElemKfi.__doc__)<jupyter_output>=====================================================
  ElemKfi(k,p,xL,xR,E,a0,a1=0,a2=0,f0=0,f1=0,f2=0)

  k :   element number (integer)
  p :   order of approximation (integer)
  xL :  x -coordinate of the left end of the element (real)
  xR :  x -coordinate of the right end of the element (real)
  E :   Young's modules ( E ) of the bar material which is taken as constant (real)
  a0,a1,a2 : coefficients in (A=a0 + a1*x + a2*(x**2)) for the area of cross-section ( A ) of the bar (real)
  f0,f1,f2 : coefficients in (f = f0 + f1*x + f2*x*x) for the distributed force ( f ) acting on the bar (real)

  Output: [Ke, fe]
  Ke : element stiffness matrix (real)
  fe : element force vector (real)
  =====================================================<jupyter_text>---  Example:1![Imgur](https://i.imgur.com/sYD9Uwt.png)<jupyter_code>[ke,fe]=ElemKfi(k=1,p=1,xL=0,xR=10,a0=6,a1=-0.1,a2=0,E=30,f0=0.1296,f1=-0.00216)
print("Local Stiffness Matrix: \n\n", ke,"\n\n")
print("Lobal Force Vector (due to distributed force):\n\n ", fe)
[ke,fe]=ElemKfi(k=1,p=1,xL=10,xR=20,a0=6,a1=-0.1,E=30,f0=0.1296,f1=-0.00216)
print("Local Stiffness Matrix: \n\n", ke,"\n\n")
print("Lobal Force Vector (due to distributed force):\n\n ", fe)
[ke,fe]=ElemKfi(k=1,p=1,xL=20,xR=30,a0=6,a1=-0.1,E=30,f0=0.1296,f1=-0.00216)
print("Local Stiffness Matrix: \n\n", ke,"\n\n")
print("Lobal Force Vector (due to distributed force):\n\n ", fe)<jupyter_output>Local Stiffness Matrix: 

 [[ 10.5 -10.5]
 [-10.5  10.5]] 


Lobal Force Vector (due to distributed force):

  [0.396 0.36 ]<jupyter_text>Example:2- book:Fundamentals of Finite Element AnalysisLinear Finite Element Analysis-Ioannis Koutromanos- example-33.2- page-87/729![Imgur](https://i.imgur.com/JlgVqxT.png) **Stiffeness and Force**<jupyter_code>[ke,fe]=ElemKfi(k=1,p=1,xL=0,xR=5,E=5000,a0=1,a1=0,a2=1/400,f0=2) #for element-1
print("Local Stiffness Matrix: \n\n", ke,"\n\n")
print("Lobal Force Vector (due to distributed force):\n\n ", fe)
[ke,fe]=ElemKfi(k=1,p=1,xL=5,xR=20,E=5000,a0=1,a1=0,a2=1/400,f0=2) #for element-2
print("Local Stiffness Matrix: \n\n", ke,"\n\n")
print("Lobal Force Vector (due to distributed force):\n\n ", fe)<jupyter_output>Local Stiffness Matrix: 

 [[ 479.16666667 -479.16666667]
 [-479.16666667  479.16666667]] 


Lobal Force Vector (due to distributed force):

  [15. 15.]<jupyter_text>Function for Global $K_{ij}$, & $f_i$ if 1-D Bar Devided in Equal Size of Elements<jupyter_code>def GlobalKfi(k,p,xL,xR,E,a0,a1=0,a2=0,f0=0,f1=0,f2=0):

  """
  =====================================================
    Note: For elements of equal length

    GlobalKfi(k,p,xL,xR,E,a0,a1=0,a2=0,f0=0,f1=0,f2=0)

    k :   number of elements (integer)
    p :   order of approximation (integer)
    xL :  x -coordinate of the left end of the element (real)
    xR :  x -coordinate of the right end of the element (real)
    E :   Young's modules ( E ) of the bar material which is taken as constant (real)
    a0,a1,a2 : coefficients in (A=a0 + a1*x + a2*(x**2)) for the area of cross-section ( A ) of the bar (real)
    f0,f1,f2 : coefficients in (f = f0 + f1*x + f2*x*x) for the distributed force ( f ) acting on the bar (real)

    Output: [Kg, fg]
    Kg : Global stiffness matrix (real)
    fg : Global force vector (real)
  =====================================================
  """

  Kg=np.zeros((k*p+1,k*p+1))   #size of element matrix
  h=xR - xL

  [n, W,xis,ShapeValueXi,dShapeValueXi]=ShapeLegRoot(p)
  fg=np.zeros(k*p+1)

  for e in range(0,k*p,p):
    #xL=N[e]
    #xR=N[e+1]
    for m in range(0,n):
      xi=xis[m]
      x = (xL+h*int(e/p))*(1 - z)/2 + (xR+h*int(e/p))*(1 + z)/2

      EA = E*(a0 + a1*x+a2*(x**2))
      EA = EA.subs(z,xi)

      f = f0 + f1*x + f2*x*x
      f=f.subs(z,xi)

      for i in range(0, p+1):
        fg[i+e] = fg[i+e] + f*ShapeValueXi[m][i]*(h/2)*W[m]

        for j in range(0,p+1):
          Kg[i+e][j+e] = Kg[i+e][j+e]+EA*dShapeValueXi[m][i]*dShapeValueXi[m][j]*(2/h)*W[m]

  return [Kg, fg]
print(GlobalKfi.__doc__)<jupyter_output>=====================================================
    Note: For elements of equal length

    GlobalKfi(k,p,xL,xR,E,a0,a1=0,a2=0,f0=0,f1=0,f2=0)

    k :   number of elements (integer)
    p :   order of approximation (integer)
    xL :  x -coordinate of the left end of the element (real)
    xR :  x -coordinate of the right end of the element (real)
    E :   Young's modules ( E ) of the bar material which is taken as constant (real)
    a0,a1,a2 : coefficients in (A=a0 + a1*x + a2*(x**2)) for the area of cross-section ( A ) of the bar (real)
    f0,f1,f2 : coefficients in (f = f0 + f1*x + f2*x*x) for the distributed force ( f ) acting on the bar (real)

    Output: [Kg, fg]
    Kg : Global stiffness matrix (real)
    fg : Global force vector (real)
  =====================================================<jupyter_text>**Example:1**![Imgur](https://i.imgur.com/sYD9Uwt.png) **Solution:**Local Stiffeness Matrices (alreday evaluated above) Global Stiffeness Matrix1. Equal Elements Length (Methode-1)<jupyter_code>[Kij, fi]=GlobalKfi(k=3,p=1,xL=0,xR=10,a0=6,a1=-0.1,a2=0,E=30,f0=0.1296,f1=-0.00216)
print("Global Stiffness Matrix: \n\n", Kij,"\n\n")
print("Global Force Vector (due to distributed force):\n\n ", fi)<jupyter_output>Global Stiffness Matrix: 

 [[ 16.5 -16.5   0.    0. ]
 [-16.5  30.  -13.5   0. ]
 [  0.  -13.5  24.  -10.5]
 [  0.    0.  -10.5  10.5]] 


Global Force Vector (due to distributed force):

  [0.612 1.08  0.864 0.36 ]<jupyter_text>---   Function for Global $K_{ij}$, & $f_i$ if 1-D Bar Devided in Unequal Size of Elements*   NodesList=Nodes X- Cordinate as list*   k=Number of element*   here,K+1=length(NodesList)<jupyter_code>def GlobalKfiU(k,p,E,NodesList,a0,a1=0,a2=0,f0=0,f1=0,f2=0):

  """
  =====================================================
    Note: For elements of unequal (or equal) length

    GlobalKfiU(k,p,E,NodesList,a0,a1=0,a2=0,f0=0,f1=0,f2=0)

    k :   number of elements (integer)
    p :   order of approximation (integer)
    NodesList:  [x0,x1,x2....xk] (X coordinates of nodes)
    E :   Young's modules ( E ) of the bar material which is taken as constant (real)
    a0,a1,a2 : coefficients in (A=a0 + a1*x + a2*(x**2)) for the area of cross-section ( A ) of the bar (real)
    f0,f1,f2 : coefficients in (f = f0 + f1*x + f2*x*x) for the distributed force ( f ) acting on the bar (real)

    Output: [Kg, fg]
    Kg : Globel stiffness matrix (real)
    fg : Globel force vector (real)
  =====================================================
  """

  Kg=np.zeros((k*p+1,k*p+1))                        #size of element matrix

  [n, W,xis,ShapeValueXi,dShapeValueXi]=ShapeLegRoot(p)
  fg=np.zeros(k*p+1)

  for e in range(0,k*p,p):
    xL=NodesList[int(e/p)]
    xR=NodesList[int(e/p)+1]
    h=xR - xL

    for m in range(0,n):
      xi=xis[m]
      x = (xL)*(1 - z)/2 + (xR)*(1 + z)/2

      EA = E*(a0 + a1*x+a2*(x**2))
      EA = EA.subs(z,xi)

      f = f0 + f1*x + f2*x*x
      f=f.subs(z,xi)

      for i in range(0, p+1):
        fg[i+e] = fg[i+e] + f*ShapeValueXi[m][i]*(h/2)*W[m]

        for j in range(0,p+1):
          Kg[i+e][j+e] = Kg[i+e][j+e]+EA*dShapeValueXi[m][i]*dShapeValueXi[m][j]*(2/h)*W[m]

  return [Kg, fg]
print(GlobalKfiU.__doc__)<jupyter_output>=====================================================
    Note: For elements of unequal (or equal) length

    GlobalKfiU(k,p,E,NodesList,a0,a1=0,a2=0,f0=0,f1=0,f2=0)

    k :   number of elements (integer)
    p :   order of approximation (integer)
    NodesList:  [x0,x1,x2....xk] (X coordinates of nodes)
    E :   Young's modules ( E ) of the bar material which is taken as constant (real)
    a0,a1,a2 : coefficients in (A=a0 + a1*x + a2*(x**2)) for the area of cross-section ( A ) of the bar (real)
    f0,f1,f2 : coefficients in (f = f0 + f1*x + f2*x*x) for the distributed force ( f ) acting on the bar (real)

    Output: [Kg, fg]
    Kg : Globel stiffness matrix (real)
    fg : Globel force vector (real)
  =====================================================<jupyter_text>**Example:2***![Imgur](https://i.imgur.com/JlgVqxT.png)* **Solution**<jupyter_code>[Kij,fi]=GlobalKfiU(k=2,p=1,NodesList=[0,5,20],E=5000,a0=1,a1=0,a2=1/400,f0=2)
print("Global Stiffness Matrix: \n\n", Kij,"\n\n")
print("Global Force Vector (due to distributed force):\n\n ", fi)<jupyter_output>Global Stiffness Matrix: 

 [[ 1020.83333333 -1020.83333333     0.        ]
 [-1020.83333333  1500.          -479.16666667]
 [    0.          -479.16666667   479.16666667]] 


Global Force Vector (due to distributed force):

  [ 5. 20. 15.]<jupyter_text>--- **Example-1 By this methode**(metode-2)<jupyter_code>[Kij,fi]=GlobalKfiU(k=3,p=1,NodesList=[0,10,20,30],a0=6,a1=-0.1,a2=0,E=30,f0=0.1296,f1=-0.00216) #global K,f
print("Global Stiffness Matrix: \n\n", Kij,"\n\n")
print("Global Force Vector (due to distributed force):\n\n ", fi)<jupyter_output>Global Stiffness Matrix: 

 [[ 16.5 -16.5   0.    0. ]
 [-16.5  30.  -13.5   0. ]
 [  0.  -13.5  24.  -10.5]
 [  0.    0.  -10.5  10.5]] 


Global Force Vector (due to distributed force):

  [0.612 1.08  0.864 0.36 ]
{'fraction_non_alphanumeric': 0.23956813539538954, 'fraction_numerical': 0.06605631747884447, 'mean_word_length': 4.339240506329114, 'pattern_counts': {'":': 0, '<': 90, '<?xml version=': 0, '>': 91, 'https://': 17, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '457907', 'license': 'MIT', 'n_tokens_mistral': 11550, 'n_tokens_neox': 10402, 'n_words': 2235, 'path': 'FEA_1D_Element.ipynb', 'repo_name': 'iampramodyadav/FEA', 'size': '27457'}