``````

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

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)))

print "  Exact Solution:\n"
print "  u = (",str(u),",",str(v),")\n"
print "  p = (",str(p),")\n"

``````
``````

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))

``````

## 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 = xsin(2piy)sin(2pix) if case == 2: u = yy(y-1) v = xx(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 [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)))

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 [ ]:

``````