In [1]:
import os, inspect
from dolfin import *
import numpy
from scipy.sparse import coo_matrix, block_diag, hstack, tril
import ExactSol
from scipy2Trilinos import scipy_csr_matrix2CrsMatrix
from PyTrilinos import Epetra, ML, AztecOO, Teuchos
import MatrixOperations as MO
import matplotlib.pylab as plt
import CheckPetsc4py as CP
import petsc4py
import sys
import HiptmairPrecond
import HiptmairApply
petsc4py.init(sys.argv)
from petsc4py import PETSc
path = os.path.abspath(os.path.join(inspect.getfile(inspect.currentframe()), ".."))
gradient_code = open(os.path.join(path, 'DiscreteGradient.cpp'), 'r').read()
compiled_gradient_module = compile_extension_module(code=gradient_code)
m = 2
errL2b =numpy.zeros((m-1,1))
errCurlb =numpy.zeros((m-1,1))
l2border = numpy.zeros((m-1,1))
Curlborder =numpy.zeros((m-1,1))
ItsSave = numpy.zeros((m-1,1))
DimSave = numpy.zeros((m-1,1))
TimeSave = numpy.zeros((m-1,1))
NN = numpy.zeros((m-1,1))
Curlgrad = numpy.zeros((m-1,1))
Massgrad = numpy.zeros((m-1,1))
Laplgrad = numpy.zeros((m-1,1))
dim = 2
for xx in xrange(1,m):
NN[xx-1] = xx+0
nn = int(2**(NN[xx-1][0]))
omega = 1
if dim == 2:
# mesh = UnitSquareMesh(int(nn),int(nn))
mesh = RectangleMesh(0.0, 0.0, 1.0, 1.5, int(nn), int(nn), 'left')
u0, p0, CurlCurl, gradPres, CurlMass = ExactSol.M2D(2,Show="yes", Mass = omega)
else:
mesh = UnitCubeMesh(int(nn),int(nn),int(nn))
u0, p0, CurlCurl, gradPres, CurlMass = ExactSol.M3D(1,Show="yes", Mass = omega)
order = 1
parameters['reorder_dofs_serial'] = False
Magnetic = FunctionSpace(mesh, "N1curl", order)
Lagrange = FunctionSpace(mesh, "CG", order)
VLagrange = VectorFunctionSpace(mesh, "CG", order)
DimSave[xx-1] = Magnetic.dim()
print Magnetic.dim()
parameters['linear_algebra_backend'] = 'uBLAS'
column = numpy.zeros(2*mesh.num_edges(), order="C") #, dtype="intc")
row = numpy.zeros(2*mesh.num_edges(), order="C") #, dtype="intc")
data = numpy.zeros(2*mesh.num_edges(), order="C") #, dtype="intc")
# Mapping = dof_to_vertex_map(Lagrange)
# c = compiled_gradient_module.Gradient(Magnetic, Mapping.astype("intc"),column,row,data)
dataX = numpy.zeros(2*mesh.num_edges(), order="C")
dataY = numpy.zeros(2*mesh.num_edges(), order="C")
dataZ = numpy.zeros(2*mesh.num_edges(), order="C")
tic()
c = compiled_gradient_module.ProlongationGradsecond(mesh, dataX,dataY,dataZ, data, row, column)
print "C++ time:", toc()
C = coo_matrix((data,(row,column)), shape=(Magnetic.dim(),Lagrange.dim())).tocsr()
G = PETSc.Mat().createAIJ(size=C.shape,csr=(C.indptr, C.indices, C.data))
Px = coo_matrix((dataX,(row,column)), shape=(Magnetic.dim(),Lagrange.dim())).tocsr()
Py = coo_matrix((dataY,(row,column)), shape=(Magnetic.dim(),Lagrange.dim())).tocsr()
Pz = coo_matrix((dataZ,(row,column)), shape=(Magnetic.dim(),Lagrange.dim())).tocsr()
Px.eliminate_zeros()
Py.eliminate_zeros()
Pz.eliminate_zeros()
if Magnetic.dim() == 8001:
VertexDoF = numpy.sin(numpy.linspace(0.0, 2*numpy.pi, num=mesh.num_vertices()))
EdgeDoFX = Px*VertexDoF
EdgeDoFY = Py*VertexDoF
EEX = Function(Magnetic)
EEY = Function(Magnetic)
VV = Function(Lagrange)
VV.vector()[:] = VertexDoF
EEX.vector()[:] = EdgeDoFX
EEY.vector()[:] = EdgeDoFY
# file = File("Plots/MagneticXdirection_"+str(Magnetic.dim())+".pvd")
# file << EEX
# file = File("Plots/MagneticYdirection_"+str(Magnetic.dim())+".pvd")
# file << EEY
# file = File("Plots/Nodal_"+str(Magnetic.dim())+".pvd")
# file << VV
plot(EEX,tite="Magnetic interpolation X-direction")
plot(EEY,tite="Magnetic interpolation Y-direction")
plot(VV,tite="Nodal represetation")
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(Magnetic,u0, boundary)
(v) = TestFunction(Magnetic)
(u) = TrialFunction(Magnetic)
(p) = TrialFunction(Lagrange)
(q) = TestFunction(Lagrange)
(Vp) = TrialFunction(VLagrange)
(Vq) = TestFunction(VLagrange)
Curl = assemble(inner(curl(u),curl(v))*dx).sparray()
Mass = assemble(inner(u,v)*dx).sparray()
Grad = assemble(inner(grad(p),v)*dx).sparray()
Laplacian = assemble(inner(grad(q),grad(p))*dx).sparray()
print "========================================"
Curlgrad[xx-1]=(Curl*C).max()
print " Curl-gradient test: ", Curlgrad[xx-1]
Massgrad[xx-1]=(Mass*C-Grad).max()
print " Mass-gradient test: ", Massgrad[xx-1]
Laplgrad[xx-1]=(Grad.transpose().tocsr()*C-Laplacian).max()
print " Lapl-gradient test: ", Laplgrad[xx-1]
print "========================================"
a = inner(curl(u),curl(v))*dx + inner(u,v)*dx
L1 = inner(v, CurlMass)*dx
Acurl,b = assemble_system(a,L1,bc)
A,b = CP.Assemble(Acurl,b)
x = b.duplicate()
ScalarLaplacian = assemble(inner(grad(p),grad(q))*dx)
VectorLaplacian = assemble(inner(grad(p),grad(q))*dx+inner(p,q)*dx)
In [10]:
e = BoundaryMesh(mesh,"interiro")
In [11]:
print e.entity_map(0)
In [ ]:
print e.entity_map
In [16]:
bc.apply(G)
In [18]:
A, N = Matrix(), 10
A.resize(N,N)
values = numpy.ones(2,dtype=numpy.float_)
rows = numpy.array([0],dtype=numpy.intc)
cols = numpy.array([1,2],dtype=numpy.intc)
A.set(values,rows,cols)
A.apply()