In [1]:

    

from dolfin import *
import numpy as np
import scipy.sparse as sp
import numpy
import matplotlib.pylab as plt
import scipy.io
from scipy2Trilinos import scipy_csr_matrix2CrsMatrix
from PyTrilinos import Epetra, ML, AztecOO, Teuchos
from dolfin import *
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
import matplotlib.pylab as plt
import PETScIO as IO
import numpy as np
import scipy.sparse as sparse
import CheckPetsc4py as CP
import scipy.sparse.linalg as sparselin
import scipy as sp
import time
nn = 2**1
mesh = RectangleMesh(0, 0, 1, 1, nn, nn,'left')

# domain_vertices = [Point(0.0, 0.0),
#                      Point(10.0, 0.0),
#                      Point(10.0, 2.0),
#                      Point(8.0, 2.0),
#                      Point(7.5, 1.0),
#                      Point(2.5, 1.0),
#                      Point(2.0, 4.0),
#                      Point(0.0, 4.0),
#                      Point(0.0, 0.0)]

#     # Create empty Mesh
# mesh = Mesh()
# PolygonalMeshGenerator.generate(mesh, domain_vertices, 0.95);
order  = 1
Magnetic = FunctionSpace(mesh, "N1curl", order)
Lagrange = FunctionSpace(mesh, "CG", order)
# L= FunctionSpace(mesh, "DG", order-1)

parameters['linear_algebra_backend'] = 'uBLAS'
b0 = Expression(("x[1]*x[1]*(x[1]-1)","x[0]*x[0]*(x[0]-1)"))


    

    




DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

In [2]:

    

C = Cell(mesh,0)


    

In [4]:

    

for cell in cells(mesh):
  print "cell", cell.index(), "has edges :", cell.entities(1)


    

    




cell 0 has edges : [0 1 2]
cell 1 has edges : [3 4 5]
cell 2 has edges : [6 7 8]
cell 3 has edges : [ 9 10 11]
cell 4 has edges : [12 13 14]
cell 5 has edges : [15 16 17]
cell 6 has edges : [18 19 20]
cell 7 has edges : [21 22 23]
cell 8 has edges : [24 25 26]
cell 9 has edges : [27 28 29]
cell 10 has edges : [30 22 31]
cell 11 has edges : [32 33 34]
cell 12 has edges : [35 36 37]
cell 13 has edges : [ 9 38 39]
cell 14 has edges : [40 41 42]
cell 15 has edges : [40  8 43]
cell 16 has edges : [44 45 46]
cell 17 has edges : [15 44 47]
cell 18 has edges : [47 48 49]
cell 19 has edges : [50 51 52]
cell 20 has edges : [12 50 53]
cell 21 has edges : [53 54 55]
cell 22 has edges : [56 57 58]
cell 23 has edges : [ 3 56 59]
cell 24 has edges : [59 60 61]
cell 25 has edges : [62 63 64]
cell 26 has edges : [18 62 65]
cell 27 has edges : [65 66 67]
cell 28 has edges : [68 69 70]
cell 29 has edges : [71 31 72]
cell 30 has edges : [73 74 75]
cell 31 has edges : [76 77 78]
cell 32 has edges : [79 80 81]
cell 33 has edges : [82 83 84]
cell 34 has edges : [85 86 87]
cell 35 has edges : [88 89 63]
cell 36 has edges : [90 91 92]
cell 37 has edges : [93 92 94]
cell 38 has edges : [95  0 96]
cell 39 has edges : [ 2 97 98]
cell 40 has edges : [ 99 100 101]
cell 41 has edges : [102  70 103]
cell 42 has edges : [104 105  87]
cell 43 has edges : [ 96  74 106]
cell 44 has edges : [107  37 106]
cell 45 has edges : [107  99 108]
cell 46 has edges : [108  97 109]
cell 47 has edges : [110 111  71]
cell 48 has edges : [ 93 112 113]
cell 49 has edges : [114 115 116]
cell 50 has edges : [117 118 119]
cell 51 has edges : [ 14 120 121]
cell 52 has edges : [ 20 122 123]
cell 53 has edges : [  5 124 125]
cell 54 has edges : [ 88 126 127]
cell 55 has edges : [ 19 127 128]
cell 56 has edges : [129 124  66]
cell 57 has edges : [129  60 122]
cell 58 has edges : [130  45 131]
cell 59 has edges : [  4 131 132]
cell 60 has edges : [133 120 134]
cell 61 has edges : [133  54 118]
cell 62 has edges : [ 17 135 136]
cell 63 has edges : [ 13 137 138]
cell 64 has edges : [139 135 137]
cell 65 has edges : [139  48  51]
cell 66 has edges : [130 140  57]
cell 67 has edges : [ 16 140 141]
cell 68 has edges : [111  77 142]
cell 69 has edges : [110 143 144]
cell 70 has edges : [ 68 145 104]
cell 71 has edges : [145 146  84]
cell 72 has edges : [147  10  75]
cell 73 has edges : [147 148 149]
cell 74 has edges : [ 35 150 151]
cell 75 has edges : [152 153 154]
cell 76 has edges : [116 155 156]
cell 77 has edges : [157 155 158]
cell 78 has edges : [159 160 113]
cell 79 has edges : [ 94 161 162]
cell 80 has edges : [ 34 163  42]
cell 81 has edges : [164 161 165]
cell 82 has edges : [159 157 164]
cell 83 has edges : [166  41 167]
cell 84 has edges : [  6 168 166]
cell 85 has edges : [169  81 170]
cell 86 has edges : [171 172 173]
cell 87 has edges : [174 175 176]
cell 88 has edges : [176 177  98]
cell 89 has edges : [169 178 152]
cell 90 has edges : [ 79 171 179]
cell 91 has edges : [180 181 182]
cell 92 has edges : [183 184 185]
cell 93 has edges : [186 187 162]
cell 94 has edges : [ 27  91 188]
cell 95 has edges : [189 186 188]
cell 96 has edges : [189 190 191]
cell 97 has edges : [192 193 149]
cell 98 has edges : [194 195 196]
cell 99 has edges : [114 195 160]
cell 100 has edges : [ 95 197 198]
cell 101 has edges : [ 73 198 103]
cell 102 has edges : [199 200  26]
cell 103 has edges : [201 200 144]
cell 104 has edges : [174 202 203]
cell 105 has edges : [204 202 205]
cell 106 has edges : [206 207 208]
cell 107 has edges : [206 209  83]
cell 108 has edges : [115 210 211]
cell 109 has edges : [ 32 212 213]
cell 110 has edges : [214 215 194]
cell 111 has edges : [214 126 210]
cell 112 has edges : [183 216 180]
cell 113 has edges : [216 217  39]
cell 114 has edges : [218 219 220]
cell 115 has edges : [218 212 221]
cell 116 has edges : [ 24 175 222]
cell 117 has edges : [223 224 146]
cell 118 has edges : [199 223 225]
cell 119 has edges : [225 226  69]
cell 120 has edges : [227 201  30]
cell 121 has edges : [227 224 228]
cell 122 has edges : [229 204  76]
cell 123 has edges : [229  25 143]
cell 124 has edges : [230 209 228]
cell 125 has edges : [208 231 134]
cell 126 has edges : [232 233   7]
cell 127 has edges : [ 90 233 234]
cell 128 has edges : [168 235 236]
cell 129 has edges : [237 235 238]
cell 130 has edges : [239 192 102]
cell 131 has edges : [239 240 105]
cell 132 has edges : [241 112 242]
cell 133 has edges : [241 234  43]
cell 134 has edges : [243   1 222]
cell 135 has edges : [243 197 226]
cell 136 has edges : [230 244  21]
cell 137 has edges : [207 244 117]
cell 138 has edges : [245 220 215]
cell 139 has edges : [245 246  89]
cell 140 has edges : [232 247 237]
cell 141 has edges : [247  28 248]
cell 142 has edges : [249 250  38]
cell 143 has edges : [249 251 148]
cell 144 has edges : [252  85  82]
cell 145 has edges : [252 253 231]
cell 146 has edges : [ 33 254 242]
cell 147 has edges : [219 254 196]
cell 148 has edges : [255 153 100]
cell 149 has edges : [255 170 151]
cell 150 has edges : [179 256 257]
cell 151 has edges : [178 256 258]
cell 152 has edges : [259  36  11]
cell 153 has edges : [259 150 260]
cell 154 has edges : [261  80 260]
cell 155 has edges : [261 172 182]

In [5]:

    

(u) = TrialFunction(Magnetic)
(v) = TestFunction(Magnetic)
(p) = TrialFunction(Lagrange)
(q) = TestFunction(Lagrange)
V = FunctionSpace(mesh, "N1curl", order)
Q = FunctionSpace(mesh, "CG", order)
W = MixedFunctionSpace([V,Q])
(uMix,pMix) = TrialFunctions(W)
(vMix,qMix) = TestFunctions(W)
a = inner(curl(v),curl(u))*dx
m = inner(u,v)*dx
b = inner(vMix,grad(pMix))*dx

# <codecell>

A = assemble(a)
M = assemble(m)
Ms = M.sparray()
A = A.sparray()

# <codecell>

B = assemble(b)
B = B.sparray()[W.dim()-V.dim():,W.dim()-Q.dim():]
ksp = PETSc.KSP().create()
# parameters['linear_algebra_backend'] = 'PETSc'
M = assemble(m)
M = CP.Assemble(M)
ksp.setOperators(M)
x = M.getVecLeft()
ksp.setFromOptions()
ksp.setType(ksp.Type.CG)
ksp.setTolerances(1e-2)
ksp.pc.setType(ksp.pc.Type.BJACOBI)

# <codecell>

OptDB = PETSc.Options()
# OptDB["pc_factor_mat_ordering_type"] = "rcm"
# OptDB["pc_factor_mat_solver_package"] = "cholmod"
ksp.setFromOptions()
C = sparse.csr_matrix((V.dim(),Q.dim()))
IO.matToSparse

# <codecell>

CC = sparse.csr_matrix((V.dim(),Q.dim()))
(v) = TrialFunction(V)
(u) = TestFunction(V)
tic()
for i in range(0,Q.dim()):
    uOut = Function(V)
    uu = Function(Q)
    x = M.getVecRight()
    zero = np.zeros((Q.dim(),1))[:,0]
    zero[i] = 1
    uu.vector()[:] = zero
    L = assemble(inner(u, grad(uu))*dx)
    rhs = IO.arrayToVec(L.array())
    ksp.solve(rhs,x)
#     x = project(grad(uu),V)
    P = x.array
    uOut.vector()[:] = P
    low_values_indices = np.abs(P) < 1e-3
    P[low_values_indices] = 0
    P=np.around(P)
    pn = P.nonzero()[0]
    for j in range(0,len(pn)):
        CC[pn[j],i] = P[pn[j]]
    del uu
print toc()


    

    




DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
/usr/local/lib/python2.7/dist-packages/scipy/sparse/compressed.py:728: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.






    




0.469259977341

In [5]:

    




    

In [2]:

    

Mmap = Magnetic.dofmap()
Lmap = Lagrange.dofmap()


    

In [3]:

    

print Lmap.vertex_to_dof_map(mesh)
print Lmap.dof_to_vertex_map(mesh)


    

    




[8 5 7 2 4 6 1 3 0]
[8 6 3 7 4 1 5 2 0]

In [28]:

    

domain_vertices = [Point(0.0, 0.0),
                     Point(10.0, 0.0),
                     Point(10.0, 2.0),
                     Point(8.0, 2.0),
                     Point(7.5, 1.0),
                     Point(2.5, 1.0),
                     Point(2.0, 4.0),
                     Point(0.0, 4.0),
                     Point(0.0, 0.0)]

    # Create empty Mesh
mesh = Mesh()
PolygonalMeshGenerator.generate(mesh, domain_vertices, 0.75);


    

In [6]:

    

edge2dof = numpy.zeros(mesh.num_edges(), dtype="int")


    

In [7]:

    

c = 0
for cell in cells(mesh):
    cellDOF = Mmap.cell_dofs(c)
    edgeVALUES = cell.entities(1)
    print edgeVALUES, cellDOF
    edge2dof[cellDOF]=edgeVALUES
    c = c+1


    

    




[0 1 2] [1 0 2]
[3 4 0] [3 4 1]
[5 4 6] [7 4 8]
[7 8 5] [11 12  7]
[ 9 10  3] [5 6 3]
[11 12  9] [ 9 10  5]
[13 12  7] [13 10 11]
[14 15 13] [14 15 13]

In [24]:

    

C = Cell(mesh,0)
edge2dof2 = edge2dof
mesh.num_edges()
C.entities(2)


    

    
Out[24]:





array([ 28,  31,  41,  74,  75, 103, 106,   8,  29,  32,  72,  73,  42,
        96,  97], dtype=uint32)

In [14]:

    

for i in range(mesh.num_cells()):
#     cell = C.end()
    C = Cell(mesh,i)
    cellDOF = Mmap.cell_dofs(i)
    edgeVALUES = C.entities(1)
    print edgeVALUES, cellDOF
    edge2dof[cellDOF]=edgeVALUES


    

    




[0 1 2] [53 47 54]
[3 4 5] [117 116 120]
[6 7 8] [172 177 173]
[ 9 10 11] [20 21 13]
[12 13 14] [89 90 86]
[15 16 17] [103 104 100]
[18 19 20] [131 132 128]
[21 22 23] [71 78 72]
[ 0 24 25] [53 56 62]
[26 27  8] [170 174 173]
[23 28 29] [72 67 73]
[30 31 32] [148 149 144]
[33 11 34] [ 9 13 14]
[ 9 35 36] [20 27 28]
[31 37 38] [149 153 154]
[39 40 41] [166 161 167]
[42 43 44] [106 109 110]
[15 42 45] [103 106 102]
[45 46 47] [102  99 105]
[48 49 50] [92 95 96]
[12 48 51] [89 92 88]
[51 52 53] [88 85 91]
[54 55 56] [114 112 115]
[ 3 54 57] [117 114 118]
[57 58 59] [118 121 122]
[60 61 62] [134 137 138]
[18 60 63] [131 134 130]
[63 64 65] [130 127 133]
[66  2 67] [61 54 60]
[66 68 28] [61 66 67]
[69 70 71] [58 63 64]
[72 73 74] [33 40 41]
[75 76 68] [75 76 66]
[33 77 78] [ 9 10  6]
[35 79 80] [27 36 37]
[81 82 21] [70 65 71]
[81 70 83] [70 63 74]
[84 61 85] [139 137 141]
[32 85 86] [144 141 145]
[87 88 89] [165 158 159]
[90 87 26] [169 165 170]
[91 74  1] [42 41 47]
[91 92 93] [42 34 43]
[94 95 96] [18 12 19]
[36 97 98] [28 25 29]
[ 99 100  69] [57 51 58]
[101 102 103] [44 38 45]
[ 72 104 105] [33 32 23]
[106 107 105] [22 15 23]
[106  94 108] [22 18 26]
[108  92 109] [26 34 35]
[110 111 112] [77 81 82]
[ 30 113 114] [148 152 147]
[115 116 117] [142 143 140]
[ 22 118 119] [78 83 84]
[ 14 120 121] [86 80 87]
[ 20 122 123] [128 126 129]
[  5 124 125] [120 123 124]
[ 84 126 117] [139 135 140]
[ 19 126 127] [132 135 136]
[128 124  64] [125 123 127]
[128  58 122] [125 121 126]
[129  43 130] [111 109 113]
[  4 130 131] [116 113 119]
[132 120  83] [79 80 74]
[132  52 118] [79 85 83]
[ 17 133 134] [100  98 101]
[ 13 135 136] [90 93 94]
[137 133 135] [97 98 93]
[137  46  49] [97 99 95]
[129 138  55] [111 107 112]
[ 16 138 139] [104 107 108]
[110 140  75] [77 68 75]
[140 141  25] [68 69 62]
[ 99 142 143] [57 52 59]
[ 82 143  67] [65 59 60]
[144  10 145] [30 21 31]
[144  79 102] [30 36 38]
[146 107  34] [11 15 14]
[146  95 147] [11 12  8]
[115 148 114] [142 146 147]
[149 150 151] [162 155 163]
[152 153  88] [157 151 158]
[ 90 154 155] [169 168 176]
[156  37  40] [160 153 161]
[156 113  89] [160 152 159]
[157 154 158] [164 168 171]
[152 149 157] [157 162 164]
[  6 159  39] [172 175 166]
[159 160 161] [175 180 181]
[162 163  77] [16 17 10]
[162 164  97] [16 24 25]
[165 166 147] [4 7 8]
[167 168 169] [1 0 2]
[167 170 165] [1 3 4]
[170 171  78] [3 5 6]
[172 173  24] [48 55 56]
[172 174  93] [48 49 43]
[175 176 155] [182 183 176]
[177 178  27] [178 179 174]
[179 175 177] [184 182 178]
[179 180 181] [184 185 186]
[182 100 101] [50 51 44]
[182 142 183] [50 52 46]
[184 150 185] [150 155 156]
[184 153 148] [150 151 146]
[186  73 183] [39 40 46]
[186 104 145] [39 32 31]

In [19]:

    

C = cells(mesh)
CC = C.next()


    

In [61]:

    

for i in range(mesh.num_edges()):
    edge = Edge(mesh,i)
normalx = numpy.zeros(mesh.num_edges())


    

In [81]:

    

for j in range(mesh.num_cells()):
    cell = Cell(mesh,j)
    edgeVALUES = cell.entities(1)
    for i in range(3):
        normalx[edgeVALUES[i]] =  cell.normal(i,0)
        print cell.normal(i,0)
        try:
            cell.normal(i,2)
            print "fail"
        except ValueError:
            print "Oops!  That was no valid number.  Try again..."
        else:
            print "2D"

print mesh.geometry()


    

    




0.707106781187
fail
2D
-1.0
fail
2D
0.0
fail
2D
0.0
fail
2D
1.0
fail
2D
-0.707106781187
fail
2D
0.707106781187
fail
2D
-1.0
fail
2D
0.0
fail
2D
0.0
fail
2D
1.0
fail
2D
-0.707106781187
fail
2D
0.707106781187
fail
2D
-1.0
fail
2D
0.0
fail
2D
0.0
fail
2D
1.0
fail
2D
-0.707106781187
fail
2D
0.707106781187
fail
2D
-1.0
fail
2D
0.0
fail
2D
0.0
fail
2D
1.0
fail
2D
-0.707106781187
fail
2D
<dolfin.cpp.mesh.MeshGeometry; proxy of <Swig Object of type 'dolfin::MeshGeometry *' at 0x40c8810> >

In [122]:

    

for i in range(mesh.num_cells()):
    cell = Cell(mesh,i)
    print cell.orientation()
#     print edge.length()


    

    




0
1
0
1
0
1
0
1

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [ ]:

    

print cell.orientation


    

In [5]:

    

(u) = TrialFunction(Magnetic)
(v) = TestFunction(Magnetic)
(p) = TrialFunction(Lagrange)
(q) = TestFunction(Lagrange)

a = inner(curl(u),curl(v))*dx 
l = inner(grad(p),grad(q))*dx+inner(p,q)*dx
# u0 = Expression(("sin(2*pi*x[1])*cos(2*pi*x[0])","-sin(2*pi*x[0])*cos(2*pi*x[1])"))
# f = 8*pow(pi,2)*u0+c*u0
CurlCurl = Expression(("-6*x[1]+2","-6*x[0]+2"))+b0
f = CurlCurl
L1  = inner(v, f)*dx


    

In [6]:

    

def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(Magnetic,b0, boundary)
Mass = assemble(inner(u,v)*dx)
bc.apply(Mass)
Acurl,b = assemble_system(a,L1,bc)
Anode = assemble(l)


    

    




DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
INFO:UFL:Adjusting missing element domain to <Domain built from <triangle cell in 2D> with label dolfin_mesh_with_id_7>.
INFO:UFL:Adjusting missing element degree to 1
INFO:UFL:Adjusting missing element domain to <Domain built from <triangle cell in 2D> with label dolfin_mesh_with_id_7>.
INFO:UFL:Adjusting missing element degree to 1
INFO:UFL:Adjusting missing element domain to <Domain built from <triangle cell in 2D> with label dolfin_mesh_with_id_7>.
INFO:UFL:Adjusting missing element degree to 1
INFO:UFL:Adjusting missing element domain to <Domain built from <triangle cell in 2D> with label dolfin_mesh_with_id_7>.
INFO:UFL:Adjusting missing element degree to 1
INFO:UFL:Adjusting missing element domain to <Domain built from <triangle cell in 2D> with label dolfin_mesh_with_id_7>.
INFO:UFL:Adjusting missing element degree to 1
INFO:UFL:Adjusting missing element domain to <Domain built from <triangle cell in 2D> with label dolfin_mesh_with_id_7>.
INFO:UFL:Adjusting missing element degree to 1
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

In [16]:

    

B = assemble(a).sparray()*C


    

    




DEBUG:FFC:Reusing form from cache.

In [17]:

    

(C-CC).todense()
# plt.show()


    

    
Out[17]:





matrix([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        ..., 
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.]])

In [64]:

    

# B.todense()


    

    
Out[64]:





matrix([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        ..., 
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.],
        [ 0.,  0.,  0., ...,  0.,  0.,  0.]])

In [ ]:

    

def build_edge2dof_map(V):
    """
    This function takes a N1Curl(1) space and return an integer valued array edge2dof.
    This array has the number of edges as its length. In particular
        edge2dof[i] = j
    means that dof #i, that is u.vector()[i], is associated to edge #j.
    """
    # Extract the cell to edge map (given an cell index, it returns the indices of its edges)
    cell2edges = V.mesh().topology()(3, 1)
    # Extract the cell dofmap (given a cell index, it returns the dof numbers) 
    cell2dofs = V.dofmap().cell_dofs
    # Array to save the result
    edge2dof = numpy.zeros(mesh.num_edges(), dtype="int")
    # Iterate over cells, associating the edges to the dofs for that cell
    for c in range(mesh.num_cells()):
        # get the global edge numbers for this cell
        c_edges = cell2edges(c)
        # get the global dof numbers for this cell
        c_dofs = cell2dofs(c)
        # associate the edge numbers to the corresponding dof numbers
        edge2dof[c_dofs] = c_edges
    # This algorithm might not look fast as it does quite some redundant work. In actual
    # runs, for most meshes, this is not the most time consuming step and does not take
    # more than a milisecond.
    return edge2dof


    

In [3]:

    

V = Magnetic
cell2edges = V.mesh().topology()(3, 1)
# Extract the cell dofmap (given a cell index, it returns the dof numbers) 
cell2dofs = V.dofmap().cell_dofs
# Array to save the result
edge2dof = numpy.zeros(mesh.num_edges(), dtype="int")
# Iterate over cells, associating the edges to the dofs fo


    

In [ ]:

    

cell2edges(1)


    

In [18]:

    

# Code for C++ evaluation of conductivity
conductivity_code = """

class Conductivity : public Expression
{
public:

  // Create expression with 3 components
  Conductivity() : Expression(3) {}

  // Function for evaluating expression on each cell
  void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
  {
    const uint D = cell.topological_dimension;
    const uint cell_index = cell.index;
    values[0] = (*c00)[cell_index];
    values[1] = (*c01)[cell_index];
    values[2] = (*c11)[cell_index];
  }

  // The data stored in mesh functions
  std::shared_ptr<MeshFunction<double> > c00;
  std::shared_ptr<MeshFunction<double> > c01;
  std::shared_ptr<MeshFunction<double> > c11;

};
"""


    

In [3]:

    

edge = Edge(mesh,0)


    

In [24]:

    

edge.entities(0)
m =edge.midpoint()
m.y()


    

    
Out[24]:





0.25

In [18]:

    

coord = mesh.coordinates()
print coord[1,:]
print coord[3,:]


    

    




[ 0.5  0. ]
[ 0.   0.5]