In [1]:

    

from dolfin import *
import numpy as np


    

In [57]:

    

parameters['reorder_dofs_serial'] = False
mesh = UnitSquareMesh(2,2)
num_cells = mesh.num_cells()
num_edges = mesh.init(1)
Q = FunctionSpace(mesh,'DG',0)
V = VectorFunctionSpace(mesh,'CG',1)
0
B = BoundaryMesh(mesh,"exterior",False)
EdgeBoundary = np.sort(B.entity_map(1).array())
print EdgeBoundary
p = Expression('x[0]')
p = interpolate(p,Q)
p = p.vector().array()
print p
v = Expression(('1','1'))
v = interpolate(v,V)
vv = v.vector().array()
print vv


    

    




DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.






    




[ 2  4  5  7 11 12 13 15]
[ 0.33333333  0.16666667  0.83333333  0.66666667  0.33333333  0.16666667
  0.83333333  0.66666667]
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]

In [59]:

    

parameters['reorder_dofs_serial'] = False
mesh = UnitSquareMesh(2,2)
num_cells = mesh.num_cells()
num_edges = mesh.init(1)
Q = FunctionSpace(mesh,'DG',0)
V = VectorFunctionSpace(mesh,'CG',1)
0
B = BoundaryMesh(mesh,"exterior",False)
EdgeBoundary = np.sort(B.entity_map(1).array())
print EdgeBoundary
p = Expression('x[0]*x[1]')
pEact = Expression('x[0]+x[1]')
p = interpolate(p,Q)
p = p.vector().array()
pEact = interpolate(pEact,Q)
pEact = pEact.vector().array()
print p
v = Expression(('1','2'))
v = interpolate(v,V)
vv = v.vector().array()
ii = 0
n = np.zeros((2,))
w = np.zeros((2,))
cellIndex = np.zeros((2,))
pJump = np.zeros((2,))
pOut = np.zeros((num_cells,))
for i in range(num_edges):
    edge = Edge(mesh,i)
#     print i, ii
    if i != EdgeBoundary[ii]:
#         print i, ii
        h = edge.length()
        v = edge.entities(0)
        v1 = Vertex(mesh,v[0])
        v2 = Vertex(mesh,v[1])
        n[0] = -(v2.x(1) - v1.x(1))
        n[1] = v2.x(0) - v1.x(0)
        n = n/np.linalg.norm(n)
        w[0] = (vv[v2.global_index()]+vv[v1.global_index()])/2
        w[1] = (vv[v2.global_index()+mesh.num_vertices()]+vv[v1.global_index()+mesh.num_vertices()])/2
        wDotn = np.dot(n,w)
        print wDotn
        j = 0
        for c in cells(edge):  
            cellIndex[j] = c.global_index()
            pJump[j] = p[c.global_index()]
            j = j+1
#         print cellIndex
        pOut[cellIndex[0]] = (1./2)*h*(pJump[0]-pJump[1])*wDotn
        pOut[cellIndex[1]] = (1./2)*h*(pJump[0]-pJump[1])*wDotn
    else:
        ii = ii + 1
#     ee = Cell(e,0)
#     print ee.global_index()
print pOut
p = Function(Q) 
p.vector()[:] = pOut
plot(p)


    

    




DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.






    




[ 2  4  5  7 11 12 13 15]
[ 0.11111111  0.02777778  0.69444444  0.44444444  0.11111111  0.02777778
  0.69444444  0.44444444]
-1.0
0.707106781187
2.0
0.707106781187
2.0
-1.0
0.707106781187
0.707106781187
[ 0.02083333 -0.04166667  0.0625     -0.125       0.02083333  0.02083333
  0.0625      0.0625    ]






    
Out[59]:





<dolfin.cpp.io.VTKPlotter; proxy of <Swig Object of type 'std::shared_ptr< dolfin::VTKPlotter > *' at 0x113e09960> >

In [42]:

    

pOut


    

    
Out[42]:





array([ 0.,  0.])

In [88]:

    

cells(1)


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-88-9decc1704a35> in <module>()
----> 1 cells(1)

/Users/michaelwathen/.hashdist/bld/profile/2vosv3gf7txo/lib/python2.7/site-packages/dolfin/cpp/mesh.py in __init__(self, *args)
    185 
    186         """
--> 187         _mesh.cells_swiginit(self,_mesh.new_cells(*args))
    188     __swig_destroy__ = _mesh.delete_cells
    189     def pos(self):

TypeError: in method 'new_cells', argument 1 of type 'dolfin::MeshEntityIteratorBase< dolfin::Cell > const &'

In [ ]:

    

e = edge.entities


    

In [ ]:

    

e = edge.mesh_id


    

In [ ]:

    

print edge.mesh


    

In [ ]:

    

print cell.cell_normal


    

In [29]:

    

mesh = UnitCubeMesh(55,1,1)
Q = FunctionSpace(mesh,'DG',0)
V = VectorFunctionSpace(mesh,'CG',1)
u_k = Function(V)
p = Expression('1')
p = interpolate(p,Q)
h = CellSize(mesh)
n = FacetNormal(mesh)
print Q.dim(), mesh.num_cells()


    

    




DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.






    




330 330

In [26]:

    

num_cells


    

    
Out[26]:





8

In [24]:

    




    

In [24]:

    




    

In [ ]: