In [26]:

    

from sympy import *
from dolfin import Expression
x = symbols('x[0]')
y = symbols('x[1]')
z = symbols('x[2]')
uu = sin(x)*exp(x+y+z)
uv = sin(y)*exp(x+y+z)
uw = sin(z)*exp(x+y+z)
# uu = (x**4+y**5+z**6)
# uv = (x**4+y**5+z**6)
# uw =(x**4+y**5+z**6)
# u = diff(uw,y)-diff(uv,z)
# v = diff(uu,z)-diff(uw,x)
# w = diff(uv,x)-diff(uu,y)
# p = sin(y)+exp(x+y+z)
u = y**1*z**2
v = x**1*z**2
w = 0
p = x*y
print "u:    ",u
print "\n"
print "v:    ",v
print "\n"
print "w:    ",w
print "\n"
print "p:    ",p
print "\n"
L1 = diff(u,x,x)+diff(u,y,y) + diff(u,z,z)
L2 = diff(v,x,x)+diff(v,y,y) + diff(v,z,z)
L3 = diff(w,x,x)+diff(w,y,y) + diff(w,z,z)
print "Laplacian 1:    ",L1
print "\n"
print "Laplacian 2:    ",L2
print "\n"
print "Laplacian 3:    ",L3
print "\n"
A1 = u*diff(u,x)+v*diff(u,y)+w*diff(u,z)
A2 = u*diff(v,x)+v*diff(v,y)+w*diff(v,z)
A3 = u*diff(w,x)+v*diff(w,y)+w*diff(w,z)

print "Advection 1:    ",A1
print "\n"
print "Advection 2:    ",A2
print "\n"
print "Advection 3:    ",A3
print "\n"
P1 = diff(p,x)
P2 = diff(p,y)
P3 = diff(p,z)
print "Pressure 1:    ",P1
print "\n"
print "Pressure 2:    ",P2
print "\n"
print "Pressure 3:    ",P3




 # (exp(x + y)*sin(y) + exp(x + y)*cos(y))**2 - 2*exp(2*x + 2*y)*sin(y)*cos(y)


    

    




u:     x[1]*x[2]**2


v:     x[0]*x[2]**2


w:     0


p:     x[0]*x[1]


Laplacian 1:     2*x[1]


Laplacian 2:     2*x[0]


Laplacian 3:     0


Advection 1:     x[0]*x[2]**4


Advection 2:     x[1]*x[2]**4


Advection 3:     0


Pressure 1:     x[1]


Pressure 2:     x[0]


Pressure 3:     0

In [18]:

    

x = symbols('x')
y = symbols('y')
z = symbols('z')


    

In [23]:

    

u.subs(z,y)


    

    
Out[23]:





x[1]*x[2]**2

In [14]:

    

u


    

    
Out[14]:





x[1]*x[2]**2

In [44]:

    

u.nsimplify()


    

    
Out[44]:





x[1]*x[2]**2

In [46]:

    

ccode(u)


    

    
Out[46]:





'x[1]*pow(x[2], 2)'

In [47]:

    

from dolfin import Expression


    

In [59]:

    

print u


    

    




x[1]*x[2]**2

In [60]:

    

x = symbols('x')
y = symbols('y')
z = symbols('z')


    

In [65]:

    

u.replace(x[0],x)


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-65-b24e2267925f> in <module>()
----> 1 u.replace(x[0],x)

TypeError: 'Symbol' object does not support indexing

In [72]:

    

u.free_symbols


    

    
Out[72]:





{x[1], x[2]}

In [ ]:

    




    

In [55]:

    

A = Matrix([u,v])


    

In [56]:

    

ccode(A)


    

    




---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-56-8404c2fde7d8> in <module>()
----> 1 ccode(A)

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/ccode.pyc in ccode(expr, assign_to, **settings)
    277 
    278     """
--> 279     return CCodePrinter(settings).doprint(expr, assign_to)
    280 
    281 

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/ccode.pyc in doprint(self, expr, assign_to)
     91                 lines.append("}")
     92         else:
---> 93             code0 = self._doprint_a_piece(expr, assign_to)
     94             lines.extend(code0)
     95 

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/codeprinter.pyc in _doprint_a_piece(self, expr, assign_to)
     37         # Setup loops over non-dummy indices  --  all terms need these
     38         if self._settings.get('contract', True):
---> 39             indices = self.get_expression_indices(expr, assign_to)
     40         else:
     41             indices = []

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/printing/codeprinter.pyc in get_expression_indices(self, expr, assign_to)
    111     def get_expression_indices(self, expr, assign_to):
    112         from sympy.tensor import get_indices, get_contraction_structure
--> 113         rinds, junk = get_indices(expr)
    114         linds, junk = get_indices(assign_to)
    115 

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/tensor/index_methods.pyc in get_indices(expr)
    235     elif expr is None:
    236         return set(), {}
--> 237     elif expr.is_Atom:
    238         return set(), {}
    239     elif isinstance(expr, Idx):

/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/matrices/matrices.pyc in __getattr__(self, attr)
   3037         else:
   3038             raise AttributeError(
-> 3039                 "%s has no attribute %s." % (self.__class__.__name__, attr))
   3040 
   3041     def integrate(self, *args):

AttributeError: MutableDenseMatrix has no attribute is_Atom.

In [31]:

    

def NS2D(case):
    x = symbols('x[0]')
    y = symbols('x[1]')

    if case == 1:
        u = sin(y)*exp(x)
        v = cos(y)*exp(x)
        p = sin(x)*cos(y)
    if case == 2:
        u = pow(y,3)
        v = pow(x,3)
        p = pow(x,2)
    if case == 3:
        uu = cos(x)*exp(x+y)
        u = diff(uu,y)
        v = -diff(uu,x)
        p = pow(x,3)*sin(y)+exp(x+y)


    L1 = diff(u,x,x)+diff(u,y,y)
    L2 = diff(v,x,x)+diff(v,y,y)

    A1 = u*diff(u,x)+v*diff(u,y)
    A2 = u*diff(v,x)+v*diff(v,y)


    P1 = diff(p,x)
    P2 = diff(p,y)


    u0 = Expression((ccode(u),ccode(v)))
    p0 = Expression(ccode(p))

    Laplacian = Expression((ccode(L1),ccode(L2)))
    Advection = Expression((ccode(A1),ccode(A2)))
    gradPres = Expression((ccode(P1),ccode(P2)))

    print "  Exact Solution:\n"
    print "  u = (",str(u),",",str(v),")\n"
    print "  p = (",str(p),")\n"
    return u0, p0, Laplacian, Advection, gradPres


    

In [ ]:

    




    

In [ ]:

    




    

In [34]:

    

ThreeD(3)


    

    




  Exact Solution:

  u = ( exp(x[0])*sin(x[1]) , exp(x[0])*cos(x[1]) )

  p = ( sin(x[0])*cos(x[1]) )







    
Out[34]:





(Coefficient(VectorElement('Lagrange', None, None, 2, None), 15),
 Coefficient(FiniteElement('Lagrange', None, None, None), 16),
 Coefficient(VectorElement('Lagrange', None, None, 2, None), 17),
 Coefficient(VectorElement('Lagrange', None, None, 2, None), 18),
 Coefficient(VectorElement('Lagrange', None, None, 2, None), 19))

In [37]:

    

def M2D(case):
    x = symbols('x[0]')
    y = symbols('x[1]')

    if case == 1:
        uu = cos(x)*exp(x+y)
        u = diff(uu,y)
        v = -diff(uu,x)
        p = x*sin(2*pi*y)*sin(2*pi*x)
    if case == 2:
        u = y*y*(y-1)
        v = x*x*(x-1)
        p = y*(y-1)*x*(x-1)

    L1 = diff(v,x,y) - diff(u,y,y)
    L2 = diff(u,x,y) - diff(v,x,x)

    P1 = diff(p,x)
    P2 = diff(p,y)

    print ccode(p)
    u0 = Expression((ccode(u),ccode(v)))
    p0 = Expression(ccode(p))

    CurlCurl = Expression((ccode(L1),ccode(L2)))
    gradPres = Expression((ccode(P1),ccode(P2)))

    print "  Exact Solution:\n"
    print "  u = (",str(u),",",str(v),")\n"
    print "  p = (",str(p),")\n"
    return u0, p0, CurlCurl, gradPres


    

In [39]:

    

p = x*sin(2*pi*y)*sin(2*pi*x)


    

In [41]:

    

s = ccode(p)


    

In [ ]: