Source code for pyOptFEM.valid2D.validStiffElas2DP1

import numpy,sympy,time
from ..FEM2D import *
from .common import *
  
def Gamma(u):
  x, y = sympy.symbols('x y')
  #return sympy.mpmath.matrix([sympy.diff(u[0],x),sympy.diff(u[1],y),sympy.diff(u[1],x)+sympy.diff(u[0],y)])
  return sympy.Matrix([sympy.diff(u[0],x),sympy.diff(u[1],y),sympy.diff(u[1],x)+sympy.diff(u[0],y)])
  

def Hooke(L,M):
  #return sympy.mpmath.matrix([[L+2*M,L,0],[L,L+2*M,0],[0,0,M]])
  return sympy.Matrix([[L+2*M,L,0],[L,L+2*M,0],[0,0,M]])
  
class TestStiffElas:
  def __init__(self, cu,cv,la,mu):
    x, y = sympy.symbols('x y')
    self.cu=cu
    self.cv=cv
    self.L=la
    self.M=mu
    self.u=sympy.sympify(cu)
    self.v=sympy.sympify(cv)
    self.fu=[sympy.Lambda((x,y),cu[0]),sympy.Lambda((x,y),cu[1])] 
    self.fv=[sympy.Lambda((x,y),cv[0]),sympy.Lambda((x,y),cv[1])] 
    if self.u[0].is_polynomial(x,y) and self.u[1].is_polynomial(x,y):
      D0=sympy.polys.Poly(self.u[0],x,y).as_dict()
      D1=sympy.polys.Poly(self.u[1],x,y).as_dict()
      self.du=max(max(sum(list(D0.keys()),1)),max(sum(list(D1.keys()),1)))
    else:
      self.du=-1
    if self.v[0].is_polynomial(x,y) and self.v[1].is_polynomial(x,y):
      D0=sympy.polys.Poly(self.v[0],x,y).as_dict()
      D1=sympy.polys.Poly(self.v[1],x,y).as_dict()
      self.dv=max(max(sum(list(D0.keys()),1)),max(sum(list(D1.keys()),1)))
    else:
      self.dv=-1
    H=Hooke(self.L,self.M)
    gU=Gamma(self.u)
    gV=Gamma(self.v)
    self.I=sympy.integrate(sympy.integrate(gV.T*H*gU,(x,0,1)),(y,0,1))[0,0]
    
def validStiffElas2DP1(**kwargs):
  """ validation of stiffness elasticity matrix assembling functions 
  
  """
  Num=kwargs.get('Num',0)
  la=kwargs.get('la',1.5)
  mu=kwargs.get('mu',0.5)
  if Num not in [0,1,2,3]:
    print("Num parameter must be in [0,1,2,3]!")
    return
      
  print('**************************************************')
  print('*     2D StiffElas Assembling P1 validations     *')
  print('**************************************************')    
  (ext,cbase)=BasesChoice(Num)
  print("  Numbering Choice : %s\n" % cbase);
  
  Th=SquareMesh(30)

  LP=[[['x - 2*y','x + y'],['x + 2*y','2*x - y'],la,mu],
      [['x**2+2*y*x+y','-2*y**2+x**2+x-y'],['3*x*y+y**2+1','3*x**2-x*y+1'],la,mu],
      [['x**3+2*y**2*x+y**2+x','y**3-2*x**2*y'],['2*x*y+y**3+x*y','3*x**3-2*x*y+x-1'],la,mu]]
      
  print('-----------------------------------------')
  print('  Test 1: Matrices errors and CPU times  ')
  print('-----------------------------------------')

  T=zeros(3)
  E=zeros(2)
  tstart=time.time()
  Mbase=StiffElasAssembling2DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.areas,la,mu,Num)
  T[0]=time.time()-tstart
  print("    Matrix size           : (%d,%d)" % Mbase.shape)
  tstart=time.time()
  MOptV1=StiffElasAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas,la,mu,Num)
  T[1]=time.time()-tstart
  E[0]=max(abs((Mbase-MOptV1).data))
  print("    Error P1base vs OptV1 : %e" % E[0])
  tstart=time.time()
  MOptV2=StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,la,mu,Num)
  T[2]=time.time()-tstart
  E[1]=max(abs((Mbase-MOptV2).data))
  print("    Error P1base vs OptV2 : %e" % E[1])
  print("    CPU times base (ref)  : %3.4f (s)" % T[0])
  print("    CPU times OptV1       : %3.4f (s) - Speed Up X%3.3f" % (T[1],T[0]/T[1]))
  print("    CPU times OptV2       : %3.4f (s) - Speed Up X%3.3f" % (T[2],T[0]/T[2]))
  checkTest1(E)

  print('-----------------------------------------------------')
  print('  Test 2: Validations by integration on [0,1]x[0,1]  ')
  print('-----------------------------------------------------')
  deg=zeros(len(LP))
  E=zeros(len(LP))
  for i in range(len(LP)):
    test=TestStiffElas(LP[i][0],LP[i][1],LP[i][2],LP[i][3])
    M=StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,test.L,test.M,Num)
    deg[i]=test.du+test.dv
    U=array([EvalFunc(test.fu[0],Th.q[:,0],Th.q[:,1]),EvalFunc(test.fu[1],Th.q[:,0],Th.q[:,1])])
    V=array([EvalFunc(test.fv[0],Th.q[:,0],Th.q[:,1]),EvalFunc(test.fv[1],Th.q[:,0],Th.q[:,1])])
    if Num==0 or Num==2:
      U=U.T;V=V.T
    U=U.reshape((2*Th.nq))
    V=V.reshape((2*Th.nq))
    Ifem=dot(M*U,V)
    E[i]=abs(Ifem-test.I)
    #print(" Ifem=%e, I=%e" % (Ifem,test.I))
    #print("Test %d: u(x,y)=%s, v(x,y)=%s" %(i,test.cu,test.cv))
    #print("  -> error : %e\n" % abs(Ifem-test.I))
    print("    function %d :\n      u(x,y)=[%s,%s],\n      v(x,y)=[%s,%s],\n           -> StiffElas error=%e" % (i,test.cu[0],test.cu[1],test.cv[0],test.cv[1],E[i]) );
  checkTest2(deg,E)

  print('--------------------------------')
  print('  Test 3: Validations by order  ')
  print('--------------------------------')
  LN=range(10,110,10)
  n=len(LN)
  h=zeros(n)
  Error=zeros(n)
  k=len(LP)-1
  test=TestStiffElas(LP[k][0],LP[k][1],LP[k][2],LP[k][3])
  print("    functions %d:\n      u(x,y)=%s,\n      v(x,y)=%s,\n      lambda=%f, mu=%f" %(k,test.cu,test.cv,test.L,test.M))
  for i in range(n):
    Th=SquareMesh(LN[i])
    meT=Th.me.T
    qT=Th.q.T;meT=Th.me.T
    tstart=time.time()
    M=StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,test.L,test.M,Num)
    TT=time.time()-tstart
    U=array([EvalFunc(test.fu[0],Th.q[:,0],Th.q[:,1]),EvalFunc(test.fu[1],Th.q[:,0],Th.q[:,1])])
    V=array([EvalFunc(test.fv[0],Th.q[:,0],Th.q[:,1]),EvalFunc(test.fv[1],Th.q[:,0],Th.q[:,1])])
    if Num==0 or Num==2:
      U=U.T;V=V.T
    U=U.reshape((2*Th.nq))
    V=V.reshape((2*Th.nq))
    Ifem=dot(M*U,V)
    h[i]=GetMaxLengthEdges(Th.q,Th.me)
    Error[i]=abs(Ifem-test.I)
    print("        Matrix size                          : (%d,%d)" % M.shape)
    print("        StiffElasAssemblingP1OptV2 CPU times : %3.3f(s)" % TT)
    print("        Error                                : %e" % Error[i]);
    #print("      Ifem=%e, I=%e" % (Ifem,test.I))
    #print("  N=%d, nq=%d, h=%.5e, Error=%.5e" % (LN[i],Th.nq,h[i],Error[i]))
    
  PlotTest3(h,Error,'Test 3 : Stiffness Elasticity Matrix (2D/P1)')